[cig-commits] [commit] knepley/feature-petsc-fe: Field: Add in PetscDS to manage discretization - Only works for a single field right now, we need to get subsetting correct - The check for the right reference cell could be improved, maybe by taking the first cell from the big DM - We are allowing the order to be determined from the command line, but we will need something else for tests (d2a8758)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/pylith

On branch  : knepley/feature-petsc-fe
Link       : https://github.com/geodynamics/pylith/compare/f33c75b19fd60eedb2a3405db76a1fee333bb1d7...5b6d812b1612809fea3bd331c4e5af98c25a536a

>---------------------------------------------------------------

commit d2a8758c9ebc398619d86e4741ab615cec795462
Author: Matthew G. Knepley <knepley at gmail.com>
Date:   Thu Oct 16 15:39:40 2014 -0500

    Field: Add in PetscDS to manage discretization
    - Only works for a single field right now, we need to get subsetting correct
    - The check for the right reference cell could be improved, maybe by taking the first cell from the big DM
    - We are allowing the order to be determined from the command line, but we will need something else for tests


>---------------------------------------------------------------

d2a8758c9ebc398619d86e4741ab615cec795462
 libsrc/pylith/topology/Field.cc | 28 ++++++++++++++++++++++++++--
 1 file changed, 26 insertions(+), 2 deletions(-)

diff --git a/libsrc/pylith/topology/Field.cc b/libsrc/pylith/topology/Field.cc
index 630afef..1ee105d 100644
--- a/libsrc/pylith/topology/Field.cc
+++ b/libsrc/pylith/topology/Field.cc
@@ -22,6 +22,8 @@
 
 #include "Mesh.hh" // USES Mesh
 
+#include "petscds.h" // USES PetscDS
+
 #include "pylith/utils/array.hh" // USES scalar_array
 
 #include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -1253,8 +1255,22 @@ pylith::topology::Field::subfieldsSetup(void)
   assert(_dm);
 
   // Setup section now that we know the total number of sub-fields and components.
-  PetscSection section = NULL;
-  PetscErrorCode err = DMGetDefaultSection(_dm, &section);PYLITH_CHECK_ERROR(err);assert(section);
+  PetscInt       dim, closureSize, vStart, vEnd, numVertices = 0;
+  PetscDS        prob      = NULL;
+  PetscFE        fe        = NULL;
+  PetscSection   section   = NULL;
+  PetscInt      *closure   = NULL, c;
+  PetscBool      isSimplex = PETSC_FALSE;
+  PetscErrorCode err;
+
+  err = DMGetDimension(_dm, &dim);PYLITH_CHECK_ERROR(err);
+  err = DMPlexGetDepthStratum(_dm, 0, &vStart, &vEnd);PYLITH_CHECK_ERROR(err);
+  err = DMPlexGetTransitiveClosure(_dm, 0, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
+  for (c = 0; c < closureSize*2; ++c) if ((closure[c] >= vStart) && (closure[c] < vEnd)) ++numVertices;
+  if (numVertices == dim+1) isSimplex = PETSC_TRUE;
+  err = DMPlexRestoreTransitiveClosure(_dm, 0, PETSC_TRUE, &closureSize, &closure);PYLITH_CHECK_ERROR(err);
+  err = DMGetDS(_dm, &prob);PYLITH_CHECK_ERROR(err);
+  err = DMGetDefaultSection(_dm, &section);PYLITH_CHECK_ERROR(err);assert(section);
   err = PetscSectionSetNumFields(section, _subfieldComps.size());PYLITH_CHECK_ERROR(err);
   err = DMSetNumFields(_dm, _subfieldComps.size());PYLITH_CHECK_ERROR(err);
   for(std::map<std::string, int>::const_iterator f_iter = _subfieldComps.begin(); f_iter != _subfieldComps.end(); ++f_iter) {
@@ -1264,7 +1280,15 @@ pylith::topology::Field::subfieldsSetup(void)
     err = PetscSectionSetFieldComponents(section, index, f_iter->second);PYLITH_CHECK_ERROR(err);
     err = DMGetField(_dm, index, &fobj);PYLITH_CHECK_ERROR(err);assert(section);
     err = PetscObjectSetName(fobj, f_iter->first.c_str());PYLITH_CHECK_ERROR(err);
+
+    std::string prefix = f_iter->first+"_";
+    /* TODO Use 0 for components of lagrange until subsetting works correctly */
+    err = PetscFECreateDefault(_dm, dim, index ? 0 : f_iter->second, isSimplex, prefix.c_str(), -1, &fe);PYLITH_CHECK_ERROR(err);
+    err = PetscObjectSetName((PetscObject) fe, f_iter->first.c_str());PYLITH_CHECK_ERROR(err);
+    err = PetscDSSetDiscretization(prob, index, (PetscObject) fe);PYLITH_CHECK_ERROR(err);
+    err = PetscFEDestroy(&fe);PYLITH_CHECK_ERROR(err);
   } // for
+  err = PetscDSSetUp(prob);PYLITH_CHECK_ERROR(err);
 
   PYLITH_METHOD_END;
 } // subfieldsSetup