MeshAdapt via Banded Interface (#788)

* Resolved logging and simmetrix compilation

* LOG: Turn on logging for the anisotropic mesh before and after

* FUNC: started adding in a hook for true model reconstruction capability

* USR: Added conditional to ensure model reconstruction only when given input flag

* ERM: Cleaned up some aspects of anisotropy.

* API: temporary fix for model reconstruction

* Anisotropy, changed to constrained volume and compute size scale and frame in same loop

* API: reconstruction function incorporated

* MAINT: cleaned out some debug statements

* ANISO: cleaned up some code for anisotropic size field

* BLD: update hashstack for updated scorec lib

* FIX: useModel incorporate into PUMIDomain

* RECONSTRUCTION: Fixed bug in loading a reconstructed mesh regarding m�material arrays by adding a necessary input flag into MeshAdaptPUMIDrvr

* RECONSTRUCTION: turn off old reconstruction call in NumericalSolution and set useModel variable if not in inputs

* RECONSTRUCTION: enable 2D reconstruction, generalization of data structures for 2D, and add default behavior to useModel flag

* RECONSTRUCTION: change the flag needed to start model reconstruction

* STACK: changed default stack hash to get updated scorec pkgs

* Update .hashstack_default

* ANISO: added maximum aspect ratio option, defaults to 100

* MA: remove input for forced adaptation

* ANISO: preserve size frame and scales post adapt

* minor typo fix in doc string

* save the meshes before and after adapt with step number

* adapt based on solve step

* Removed clamping from anisotropy, gradeMesh instead of smoothing

* Changed how the adapt inputs are passed in to preserve input fields

* Added feature to search for error near interface when no error target is supplied, turned clamping back on

* anisotropic grading added through aspect ratio scaling.

* Use the minimum height function and change the error to size relationship

* error-to-size relationship fix for isotropic, setup framework for anisotropic logic

* FIX: fix some errors with API

* FIX: resolve test for gmshCheck by adding in targetError

* API: added several inputs/parameters to MeshAdaptPUMIDrvr including hPhi, numAdaptSteps (for controlling adaptation) in NumericalSolution

* FIX: avoid double counting of solve step when counting for adaptivity

* API: enable transfer of epsFact into size field considerations

* STATE: save state of initial implementation of parallel gradation algorithm

* FEATURE: implementation of parallel gradation

* FIX/API: require the use of context in order to get domain and physics

* FIX: fix merge conflict

* FIX: fix conflicts from merge

* FIX:avoid node material conflicts by creating a field beforehand and using that as the source of classification

* MAINT: modified some pList references to use context, ct // turned off goodQuality input

* FIX: resolve parallel assignment of node materials through parallel communication

* FIX/MAINT: cleanup of debug statements, added description, simplify nodeMaterial procedure to resolve bug

* OPT: move mesh derivation routine to immediately after mesh generation

* ENH: improve reliability of updateMaterial from a derived model by taking the minimum material tag

* MAINT: for more consistent behavior on HPC, use Zoltan only for now

* MAINT: turn off mass check since not all applications have the requisite fields

* FIX: correct numbering of logging meshes to prevent overlapping names in parallel cases

* MAINT: removed some old code/clutter

* ENH: Added edge-based band intersection detection for interface based refinement, accounting for parallel cases as well

* MAINT/DOC: cleaned up debug statements

* MAINT/ENH: Added initial adapt before solve and moved some code into a function

* MAINT: turn adapt back on, remove unnecessary statements after adapt in initial adapt

* SAVE: state of numericalsolution with list of needed data structures to transfer/rebuild

* FIX: ensure consistency in adapt workflow by setting free DOFs from RD and MCorr

* FIX: matched shock capturing for VOF by transporting old dofs and performing computations

* Save state for debugging, single phase vbdf time integration with adaptivity workflow is consistent

* FIX: generalized postAdapt functions as member functions of relevant classes,can handle multiphase 2D vbdf case with systemStepExact=True.

* Turned off vof shock capturing lagging. Need to store u^{n-2} solution to get it working properly

* ENH/FIX: Added secondary initial adapt and corrected how nUpdatesTimeHistoryCalled is updated after adapt

* FIX: avoid memory error by properly destorying marker field

* TRY/FIX: force shock capturing lagging post adapt based on incorrect numDiff last and rectify subgrid error array for lagging

* UPDATE: hard-coded parameter for interface band size

* FIX: wrap initial adapt in sfconfig dependent conditional

* TESTS: fix tests that require epsFact information to be passed in

* TEST: didn't finish fixing test from previous commit

* FIX: needed to destroy some entities to clear memory

* VMS: finish parallel VMS implementation with simple all_reduce sum. Also extension of centroid computation to 3D.

* FIX: Add definition of postAdaptUpdate for backward euler cfl

* MAINT: add warning if users ever call postAdaptUpdate() without having a proper definition

* FIX: shock capturing depends on subgrid error so sge needs to be in loop when recovering the shock capturing

* FIX: resolved issues with lagging with multiple models by defining appropriate stages that represent the model stagger. Also ensured adaptive time stepping is accounted for with the dtLast modification.

* MAINT: consistency in function def

* FIX: MCorr no longer needs to be corrected since LS ebqe array was fixed

* FIX: account for situation when dtLast does not exist

* FIX: is not instead of !=

* FIX: populate the lagged eddy-viscosity array. This for some reason is not sufficient to yield exactly the same results between pseudoadapt and non-adaptive case. 1e-12 error.

* TST: added pseudo option for adapt to enable ease of testing

* TST: added gauge comparison test for post-adapt restart/solution transfer via 2D dambreak case

* TST: added gauge compare test to setup.py

* TEST: hack to get really crude test of gauge comparison running. will need to reformat into template of other tests

* FIX/TST: gauge compare test not being run by Travis

* FIX/TST: change the install path for gauge_compare test to fix Travis error

* DOC: add some documentation strings to the grading functions

* DOC: added some more documentation strings

.gitignore

@@ -22,6 +22,7 @@ data*txt
 *.poly
 *.bak
 *.log.*
+*.swp
 install_*
 config_*
 build/

proteus/MeshAdaptPUMI/MeshAdaptPUMI.h

@@ -18,7 +18,7 @@
 class MeshAdaptPUMIDrvr{
 
   public:
-  MeshAdaptPUMIDrvr(double, double, int, const char*, const char*,const char*,double,double,int,double,double); 
+  MeshAdaptPUMIDrvr(double, double, double, int, int, int, const char*, const char*,const char*,double,double,int,double,double); 
   ~MeshAdaptPUMIDrvr();
 
   int loadModelAndMesh(const char* modelFile, const char* meshFile); //load the model and mesh
@@ -47,7 +47,7 @@ class MeshAdaptPUMIDrvr{
   //Functions used to transfer information between PUMI and proteus
   int transferFieldToPUMI(const char* name, double const* inArray, int nVar, int nN);
   int transferFieldToProteus(const char* name, double* outArray, int nVar, int nN);
-  int transferPropertiesToPUMI(double* rho_p, double* nu_p,double* g_p, double deltaT);
+  int transferPropertiesToPUMI(double* rho_p, double* nu_p,double* g_p, double deltaT, double interfaceBandSize);
   //int transferBCtagsToProteus(int* tagArray, int idx, int* ebN, int* eN_global, double* fluxBC);
   //int transferBCsToProteus();
 
@@ -70,12 +70,15 @@ class MeshAdaptPUMIDrvr{
 
 
   //Public Variables
-  double hmax, hmin; //bounds on mesh size
+  double hmax, hmin, hPhi; //bounds on mesh size
   int numIter; //number of iterations for MeshAdapt
   int nAdapt; //counter for number of adapt steps
   int nEstimate; //counter for number of error estimator calls
   int nsd; //number of spatial dimensions
   int maxAspect; //maximum aspect ratio
+  int adaptMesh; //is adaptivity turned on?
+  int numAdaptSteps; //Number adaptivity
+  double N_interface_band; //number of elements in half-band around interface
   double gradingFactor;
 
   //User Inputs

proteus/MeshAdaptPUMI/MeshAdaptPUMI.pyx

@@ -17,9 +17,10 @@ cdef extern from "cmeshToolsModule.h":
 
 cdef extern from "MeshAdaptPUMI/MeshAdaptPUMI.h":
     cdef cppclass MeshAdaptPUMIDrvr:
-        MeshAdaptPUMIDrvr(double, double, int, char*, char*,char*,double,double,int,double,double)
+        MeshAdaptPUMIDrvr(double, double, double, int, int, int, char*, char*,char*,double,double,int,double,double)
         int numIter, numAdaptSteps
         string size_field_config, adapt_type_config
+        int adaptMesh
         int isReconstructed
         int loadModelAndMesh(char *, char*)
         int getSimmetrixBC()
@@ -32,7 +33,7 @@ cdef extern from "MeshAdaptPUMI/MeshAdaptPUMI.h":
         int updateMaterialArrays2(Mesh&)
         int transferFieldToPUMI(char*, double*, int, int)
         int transferFieldToProteus(char*, double*, int, int)
-        int transferPropertiesToPUMI(double*, double*,double*,double)
+        int transferPropertiesToPUMI(double*, double*,double*,double,double)
         int transferModelInfo(int*,int*,int*,int*,int*,int*,int)
         int transferBCtagsToProteus(int*, int, int*, int*,double*)
         int transferBCsToProteus()
@@ -49,13 +50,13 @@ cdef extern from "MeshAdaptPUMI/MeshAdaptPUMI.h":
 
 cdef class MeshAdaptPUMI:
     cdef MeshAdaptPUMIDrvr *thisptr
-    cdef double hmax, hmin
+    cdef double hmax, hmin, hPhi
     cdef int numIter, numAdaptSteps
     cdef int isReconstructed
-    cdef double deltaT
-    def __cinit__(self, hmax=100.0, hmin=1e-8, numIter=10, sfConfig="ERM",maType="test",logType="off",targetError=0,targetElementCount=0,reconstructedFlag=0,maxAspectRatio=100.0,gradingFact=1.5):
+    cdef int adaptMesh 
+    def __cinit__(self, hmax=100.0, hmin=1e-8, hPhi=1e-2, adaptMesh=0, numIter=10, numAdaptSteps=10 ,sfConfig="ERM",maType="test",logType="off",targetError=0,targetElementCount=0,reconstructedFlag=0,maxAspectRatio=100.0,gradingFact=1.5):
         logEvent("MeshAdaptPUMI: hmax = {0} hmin = {1} numIter = {2}".format(hmax,hmin,numIter))
-        self.thisptr = new MeshAdaptPUMIDrvr(hmax, hmin, numIter, sfConfig,maType,logType,targetError,targetElementCount,reconstructedFlag,maxAspectRatio,gradingFact)
+        self.thisptr = new MeshAdaptPUMIDrvr(hmax, hmin, hPhi, adaptMesh, numIter, numAdaptSteps, sfConfig,maType,logType,targetError,targetElementCount,reconstructedFlag,maxAspectRatio,gradingFact)
     def __dealloc__(self):
         del self.thisptr
     def size_field_config(self):
@@ -64,6 +65,10 @@ cdef class MeshAdaptPUMI:
         return self.thisptr.adapt_type_config
     def numIter(self):
         return self.thisptr.numIter
+    def adaptMesh(self):
+        return self.thisptr.adaptMesh
+    def numAdaptSteps(self):
+        return self.thisptr.numAdaptSteps
     def isReconstructed(self):
         return self.thisptr.isReconstructed
     def loadModelAndMesh(self, geomName, meshName):
@@ -99,11 +104,11 @@ cdef class MeshAdaptPUMI:
     def transferFieldToProteus(self, name, np.ndarray[np.double_t,ndim=2,mode="c"] outArray):
         outArray = np.ascontiguousarray(outArray)
         return self.thisptr.transferFieldToProteus(name, &outArray[0,0], outArray.shape[1], outArray.shape[0])
-    def transferPropertiesToPUMI(self, np.ndarray[np.double_t,ndim=1,mode="c"] rho, np.ndarray[np.double_t,ndim=1,mode="c"] nu, np.ndarray[np.double_t,ndim=1,mode="c"] g, deltaT):
+    def transferPropertiesToPUMI(self, np.ndarray[np.double_t,ndim=1,mode="c"] rho, np.ndarray[np.double_t,ndim=1,mode="c"] nu, np.ndarray[np.double_t,ndim=1,mode="c"] g, double deltaT, double interfaceBandSize):
         rho = np.ascontiguousarray(rho)
         nu = np.ascontiguousarray(nu)
         g = np.ascontiguousarray(g)
-        return self.thisptr.transferPropertiesToPUMI(&rho[0],&nu[0],&g[0],deltaT)
+        return self.thisptr.transferPropertiesToPUMI(&rho[0],&nu[0],&g[0],deltaT,interfaceBandSize)
     def transferModelInfo(self, np.ndarray[int,ndim=1,mode="c"] numModelEntities,
                                 np.ndarray[int,ndim=2,mode="c"] edges,
                                 np.ndarray[int,ndim=2,mode="c"] faces,

proteus/MeshAdaptPUMI/MeshConverter.cpp

@@ -1,6 +1,7 @@
 #include <algorithm>
 
 #include "MeshAdaptPUMI.h"
+#include <PCU.h>
 #include "mesh.h"
 #include <apfShape.h>
 
@@ -513,26 +514,98 @@ int MeshAdaptPUMIDrvr::updateMaterialArrays2(Mesh& mesh)
   apf::MeshIterator* it;
   apf::MeshEntity* f;
 
+  //first associate all nodes with a material tag and synchronize fields to avoid mismatches 
+  //The procedure is to have each vertex look for its classification.
+  //If it is classified in the region, then it is interior.
+  //Else, loop over adjacent faces and stop at first instance of mesh face classified on model boundary and take tag.
+    //If there are no such adjacent mesh faces, then set the value to be -1. This should only happen if the vertex is a shared entity.
+  //If the vertex is shared, communicate value to remote copies.
+  //When receiving values, if the current value is -1, write to field the received value. Otherwise, do nothing. 
+
+  apf::Field* nodeMaterials = apf::createLagrangeField(m, "nodeMaterials", apf::SCALAR, 1);
+  it = m->begin(0);
+  PCU_Comm_Begin();
+  while(f = m->iterate(it))
+  {
+    geomEnt = m->toModel(f);
+    //if classified in a region
+    if(m->getModelType(geomEnt) == m->getDimension())
+    {
+      apf::setScalar(nodeMaterials,f,0,0); 
+    }
+    else
+    {
+      apf::Adjacent vert_adjFace;
+      m->getAdjacent(f,m->getDimension()-1,vert_adjFace);
+      apf::MeshEntity* face;
+      for(int i =0; i<vert_adjFace.getSize();i++)
+      {
+        face=vert_adjFace[i];
+        geomEnt = m->toModel(face);
+
+        //IF mesh face is classified on boundary
+        if(m->getModelType(geomEnt) == m->getDimension()-1)
+        {
+          geomTag = m->getModelTag(geomEnt);
+          apf::setScalar(nodeMaterials,f,0,geomTag);
+          if(m->isShared(f))
+          {
+            apf::Copies remotes;
+            m->getRemotes(f,remotes);
+            for(apf::Copies::iterator iter = remotes.begin(); iter != remotes.end(); ++iter)
+            {
+              PCU_COMM_PACK(iter->first,iter->second);
+              PCU_COMM_PACK(iter->first,geomTag);
+            }
+          }  
+          break;
+        }
+        if(i == vert_adjFace.getSize()-1 )
+          apf::setScalar(nodeMaterials,f,0,-1);
+      }
+    }
+  }
+  PCU_Comm_Send();
+  while(PCU_Comm_Receive())
+  {
+    PCU_COMM_UNPACK(f);
+    PCU_COMM_UNPACK(geomTag);
+    int currentTag = apf::getScalar(nodeMaterials,f,0);
+    int newTag = std::min(currentTag,geomTag);
+    //if vertex is not interior and had no adjacent faces, take received values
+    //else take minimum value of all tags
+    if(currentTag==-1)
+      apf::setScalar(nodeMaterials,f,0,geomTag);
+    else
+      apf::setScalar(nodeMaterials,f,0,newTag);
+  }
+  //Ensure there are no mismatches across parts and then assign node materials
+  apf::synchronize(nodeMaterials);
+  it = m->begin(0);
+  while(f=m->iterate(it))
+  {
+    int vID = localNumber(f);
+    mesh.nodeMaterialTypes[vID] = apf::getScalar(nodeMaterials,f,0);
+  }
+
   //First iterate over all faces in 3D, get the model tag and apply to all downward adjacencies
   int dim = m->getDimension()-1;
   it = m->begin(dim);
-  while(f = m->iterate(it)){
+  while(f = m->iterate(it))
+  {
     int i = localNumber(f);
     geomEnt = m->toModel(f);
     geomTag = m->getModelTag(geomEnt);
-    if(m->getModelType(geomEnt) == dim){
+    if(m->getModelType(geomEnt) == dim)
+    {
       mesh.elementBoundaryMaterialTypes[i] = geomTag;
-      apf::Adjacent face_adjVert;
-      m->getAdjacent(f,0,face_adjVert);
-      for(int j=0;j<face_adjVert.getSize();j++){
-          int vID = localNumber(face_adjVert[j]);
-          mesh.nodeMaterialTypes[vID] = geomTag;
-      }
     }
   }
   m->end(it);
 
+  apf::destroyField(nodeMaterials);
   std::cout<<"Finished faces and verts\n";
+
   //Loop over regions
   dim = m->getDimension();
   it = m->begin(dim);
@@ -555,7 +628,7 @@ int MeshAdaptPUMIDrvr::updateMaterialArrays2(Mesh& mesh)
  * scorec/core. This may be added into scorec/core eventually and removed.
  */
 
-#include <PCU.h>
+//#include <PCU.h>
 #include "apfConvert.h"
 #include "apfMesh2.h"
 #include "apf.h"

proteus/MeshAdaptPUMI/MeshFields.cpp

@@ -114,7 +114,7 @@ int MeshAdaptPUMIDrvr::transferFieldToProteus(const char* name, double* outArray
   return 0;
 }
 
-int MeshAdaptPUMIDrvr::transferPropertiesToPUMI(double* rho_p, double* nu_p, double *g_p, double deltaT)
+int MeshAdaptPUMIDrvr::transferPropertiesToPUMI(double* rho_p, double* nu_p, double *g_p, double deltaT, double interfaceBandSize)
 /**
  * @brief Transfer material properties to the MeshAdaptPUMI class
  * 
@@ -134,7 +134,8 @@ int MeshAdaptPUMIDrvr::transferPropertiesToPUMI(double* rho_p, double* nu_p, dou
   }
   else if(nsd==3){
     g[0] = g_p[0]; g[1] = g_p[1]; g[2] = g_p[2];
-  } 
+  }
+  N_interface_band = interfaceBandSize;
   delta_T = deltaT;
   return 0;
 }