Merging develop into master

Author: Nightly Bot

Diff for: apf_sim/CMakeLists.txt

@@ -27,6 +27,8 @@ set(USE_FIELDSIM FALSE)
 if( ${SIMMODSUITE_SimField_FOUND} AND ENABLE_FIELDSIM )
   set(USE_FIELDSIM TRUE)
 endif()
+set(USE_SIM_ADVMESHING ${HAVE_SIMADVMESHING})
+
 configure_file(
     "${CMAKE_CURRENT_SOURCE_DIR}/apf_simConfig.h.in"
     "${CMAKE_CURRENT_BINARY_DIR}/apf_simConfig.h")

Diff for: apf_sim/apfSIM.cc

@@ -896,6 +896,8 @@ void MeshSIM::writeNative(const char* fileName)
 
 void MeshSIM::destroyNative()
 {
+  while(this->countFields())
+    apf::destroyField(this->getField(0));
   M_release(mesh);
 }
 

Diff for: apf_sim/apfSIMDataOf.h

@@ -20,6 +20,11 @@ class SIMDataOf : public FieldDataOf<T>
     fd = fd_in;
     pf = static_cast<pPolyField>(Field_def(fd));
   }
+  ~SIMDataOf()
+  {
+    Field_release(fd);
+    FieldDef_delete(pf);
+  }
   virtual void init(FieldBase * f)
   {
     FieldData::field = f;

Diff for: apf_sim/apf_simConfig.h.in

@@ -1 +1,2 @@
 #cmakedefine USE_FIELDSIM
+#cmakedefine USE_SIM_ADVMESHING

Diff for: gmi_sim/gmi_sim.cc

@@ -193,20 +193,22 @@ static gmi_set* edge_regions(pGEdge e)
   pPList list = GE_faces((pGEdge)e);
 
   std::set<pGRegion> rgns;
-  gmi_set* regions_set;
   for (int i = 0; i < PList_size(list); i++) {
-   regions_set = face_regions((pGFace)PList_item(list, i));
+   gmi_set* regions_set = face_regions((pGFace)PList_item(list, i));
    for (int j = 0; j < regions_set->n; j++) {
      rgns.insert( (pGRegion) regions_set->e[j] );
    }
+   gmi_free_set(regions_set);
   }
+  PList_delete(list);
 
   int count = 0;
   gmi_set* s = gmi_make_set(rgns.size());
   for (std::set<pGRegion>::iterator it=rgns.begin(); it!=rgns.end(); ++it) {
     s->e[count] = (gmi_ent*)*it;
     count++;
   }
+  rgns.clear();
 
   return s;
 }

Diff for: phasta/CMakeLists.txt

@@ -43,6 +43,7 @@ if(ENABLE_SIMMETRIX)
       phAttrib.cc
       phSnap.cc
       phMeshQuality.cc
+      phRigidBody.cc
      )
 endif()
 
@@ -62,6 +63,8 @@ set(HEADERS
   phstream.h
   phiotimer.h
   phastaChef.h
+  phBC.h
+  ph.h
 )
 
 # Add the ph library

Diff for: phasta/phGeomBC.cc

@@ -335,6 +335,7 @@ void writeGeomBC(Output& o, std::string path, int timestep)
     params[1] = 3;
     ph_write_ints(f, "m2g classification", o.arrays.m2gClsfcn,
       params[0] * params[1], 2, params);
+    params[0] = m->count(0);
     params[1] = 2;
     ph_write_doubles(f, "m2g parametric coordinate", o.arrays.m2gParCoord,
       params[0] * params[1], 2, params);

Diff for: phasta/phGrowthCurves.cc

@@ -73,6 +73,7 @@ void getGrowthCurves(Output& o)
           break;
         }
       }
+      VIter_delete(vIter);
 
       if(isBoundaryLayerFace){
 //      Add gFace to gEntities
@@ -104,6 +105,7 @@ void getGrowthCurves(Output& o)
         }
       }
     }
+    GFIter_delete(gFIter);
 
 //  get seeds of all growth curves
     pPList allSeeds = PList_new();
@@ -175,6 +177,7 @@ void getGrowthCurves(Output& o)
         //Append seeds to allSeeds
         PList_appPList(allSeeds,seeds);
       }
+      VIter_delete(vIter);
     }
 
 //  get info of growth curves

Diff for: phasta/phInput.cc

@@ -79,6 +79,8 @@ static void setDefaults(Input& in)
   in.betaSize = 2e-5;
   in.gammaDist = 3e-6;
   in.gammaSize = 3e-5;
+  in.nRigidBody = 0;
+  in.nRBParam = 12;
 }
 
 Input::Input()
@@ -157,6 +159,8 @@ static void formMaps(Input& in, StringMap& stringMap, IntMap& intMap, DblMap& db
   dblMap["betaSize"] = &in.betaSize;
   dblMap["gammaDist"] = &in.gammaDist;
   dblMap["gammaSize"] = &in.gammaSize;
+  intMap["nRigidBody"] = &in.nRigidBody;
+  intMap["nRBParam"] = &in.nRBParam;
 }
 
 template <class T>

Diff for: phasta/phInput.h

@@ -11,6 +11,7 @@
 */
 
 #include <string>
+#include <vector>
 
 struct RStream;
 
@@ -170,6 +171,12 @@ class Input
     double gammaDist;
     /** \brief absolute isotropic size within [betaDist:gammDist) of zero level set */
     double gammaSize;
+    /** \brief number of rigid bodies */
+    int nRigidBody;
+    /** \brief number of parameters for each rigid body */
+    int nRBParam;
+    /** \brief parameter data for rigid body */
+    std::vector<double> rbParamData;
 };
 
 int countNaturalBCs(Input& in);

Diff for: phasta/phMeshQuality.cc

@@ -8,40 +8,61 @@
 #include <apfMatrix.h>
 #include <maSize.h>
 #include <maShape.h>
+//#include "apf_simConfig.h"
+
+#ifdef HAVE_SIMMETRIX
 #include <SimUtil.h>
 #include <SimModel.h>
 #include <SimPartitionedMesh.h>
 #include <SimMeshTools.h>
+//#ifdef HAVE_SIMADVMESHING
+#include <SimAdvMeshing.h>
+//#endif
+#endif
 
+#include <PCU.h>
 #include <pcu_util.h>
 #include <cstdio>
 
 namespace ph {
 
+static void initArray(double array[3][3]) {
+  for(int i = 0; i < 3; i++) {
+    for(int j = 0; j < 3; j++) {
+      array[i][j] = 0.0;
+    }
+  }
+}
+
 void attachSIMSizeField(apf::Mesh2* m, apf::Field* sf_mag, apf::Field* sf_dir) {
 // loop over all vertices
   double size[1];
   double anisosize[3][3];
   apf::MeshEntity* v;
   apf::MeshIterator* vit = m->begin(0);
   while ((v = m->iterate(vit))) {
+#ifdef HAVE_SIMMETRIX
 // get sim size field
 // this is for simmetrix mesh, should be generalized
-    pVertex meshVertex = reinterpret_cast<pVertex>(v);
-    int sztype = V_size(meshVertex, size, anisosize);
-    PCU_ALWAYS_ASSERT(sztype == 1 || sztype == 2);
-
-// consider iso size field as anisotropic size field
-    if (sztype == 1) {
-      for(int i = 0; i < 3; i++) {
-        for(int j = 0; j < 3; j++) {
-          if (i==j)
-            anisosize[i][j] = size[0];
-          else
-            anisosize[i][j] = 0.0;
-        }
+    {
+      pVertex meshVertex = reinterpret_cast<pVertex>(v);
+      int sztype = V_size(meshVertex, size, anisosize);
+      PCU_ALWAYS_ASSERT(sztype == 1 || sztype == 2);
+
+      // consider iso size field as anisotropic size field
+      if (sztype == 1) {
+        initArray(anisosize);
+        anisosize[0][0] = size[0];
+        anisosize[1][1] = size[0];
+        anisosize[2][2] = size[0];
       }
     }
+#else
+    {
+      PCU_ALWAYS_ASSERT_VERBOSE(0,
+         "please turn on Simmetrix and re-compile the code.\n");
+    }
+#endif
 
 // transfer to apf field sf_mag and sf_dir
 /* note that the frame in Simmetrix is stored by row
@@ -76,14 +97,22 @@ extern"C"{
 /* input x1,x2,x3 are current coordinates of mesh in phasta
    output minq is the minimum quality of mesh */
 void core_measure_mesh (double x1[], double x2[], double x3[], int numnp,
-                        double& minq) {
+                        double& minvq, double& minfq) {
+  double** op = (double**)malloc(sizeof(double*) * numnp);
+  for (int i = 0; i < numnp; i++)
+    op[i] = (double*)malloc(sizeof(double) * 3);
 // loop over all vertices
   int counter = 0;
+  apf::Vector3 p;
   apf::MeshEntity* v;
   apf::MeshIterator* vit = m->begin(0);
   while ((v = m->iterate(vit))) {
+// get original mesh coordinates
+    m->getPoint(v, 0, p);
+    op[counter][0] = p[0];
+    op[counter][1] = p[1];
+    op[counter][2] = p[2];
 // update the coordinates of current mesh
-    apf::Vector3 p;
     p[0] = x1[counter];
     p[1] = x2[counter];
     p[2] = x3[counter];
@@ -96,24 +125,107 @@ void core_measure_mesh (double x1[], double x2[], double x3[], int numnp,
 // transfer to apf sizefield
   if(m->findField("sizes")) apf::destroyField(m->findField("sizes"));
   if(m->findField("frames")) apf::destroyField(m->findField("frames"));
-// this is for simmetrix mesh, should be generalized
   apf::Field* sf_mag = apf::createSIMFieldOn(m, "sizes", apf::VECTOR);
   apf::Field* sf_dir = apf::createSIMFieldOn(m, "frames", apf::MATRIX);
+// this is for simmetrix mesh, should be generalized
+  PCU_ALWAYS_ASSERT_VERBOSE(in.simmetrixMesh, "core_measure_mesh only supports Simmetrix mesh!\n");
   ph::attachSIMSizeField(m, sf_mag, sf_dir);
   ma::SizeField* sf = ma::makeSizeField(m, sf_mag, sf_dir);
 
-// measure
-  minq = 1.0;
+// measure mesh region
+  minvq = 1.0;
   apf::MeshEntity* e;
   apf::MeshIterator* eit = m->begin(m->getDimension());
   while( (e = m->iterate(eit)) ) {
     if (! m->isOwned(e))
       continue;
-    double q = ma::measureElementQuality(m, sf, e);
-    if (q < minq)
-      minq = q;
+    /* not support other mesh region type */
+    if (m->getType(e) != apf::Mesh::TET)
+      continue;
+#ifdef HAVE_SIMMETRIX
+    if (EN_isBLEntity(reinterpret_cast<pEntity>(e)))
+      continue;
+#endif
+    double vq = ma::measureElementQuality(m, sf, e); // cubic mean ratio
+    vq = cbrt(vq); // mean ratio
+    if (vq < minvq)
+      minvq = vq;
   }
   m->end(eit);
+
+// measure mesh face in layered mesh
+  minfq = 1.0;
+#ifdef HAVE_SIMMETRIX
+  if (in.simmetrixMesh == 1) {
+// get simmetrix mesh
+    apf::MeshSIM* apf_msim = dynamic_cast<apf::MeshSIM*>(m);
+    pParMesh pmesh = apf_msim->getMesh();
+    pMesh mesh = PM_mesh(pmesh,0);
+
+// get simmetrix model
+    gmi_model* gmiModel = apf_msim->getModel();
+    pGModel model = gmi_export_sim(gmiModel);
+
+// measure triangles
+    pGFace gFace;
+    pFace meshFace;
+    pEntity seedRegion;
+    pPList growthRegion = PList_new();
+    pPList growthLayerFace = PList_new();
+    GFIter gFIter = GM_faceIter(model);
+    while((gFace = GFIter_next(gFIter))){
+      FIter fIter = M_classifiedFaceIter(mesh, gFace, 1);
+      while((meshFace = FIter_next(fIter))){
+        if(BL_isBaseEntity(meshFace,gFace) == 1){
+          for (int fromSide = 0; fromSide < 2; fromSide++) {
+            int hasSeed = BL_stackSeedEntity(meshFace,gFace,fromSide,NULL,&seedRegion);
+            PCU_ALWAYS_ASSERT_VERBOSE(hasSeed >= 0, "BL blending is not supported!\n");
+            if (hasSeed == 0)
+              continue;
+            PCU_ALWAYS_ASSERT(BL_growthRegionsAndLayerFaces
+                      ((pRegion)seedRegion,growthRegion,growthLayerFace,Layer_Entity) == 1);
+            for (int iglf = 0; iglf < PList_size(growthLayerFace); iglf++) {
+              pFace blFace = (pFace)PList_item(growthLayerFace, iglf);
+              double fq = ma::measureElementQuality(m, sf, reinterpret_cast<apf::MeshEntity*> (blFace)); // squared mean ratio
+              fq = (fq > 0) ? sqrt(fq) : -sqrt(-fq); // mean ratio
+              if (fq < minfq)
+                minfq = fq;
+            }
+            PList_clear(growthRegion);
+            PList_clear(growthLayerFace);
+          }
+        }
+      }
+      FIter_delete(fIter);
+    }
+    GFIter_delete(gFIter);
+    PList_delete(growthRegion);
+    PList_delete(growthLayerFace);
+  }
+  else
+#endif
+  {
+    printf("PUMI-based mesh doesn't have BL info. minfq is set to be 1.0\n");
+  }
+
+// restore mesh
+  counter = 0;
+  vit = m->begin(0);
+  while ((v = m->iterate(vit))) {
+// use the original coordinates
+    p[0] = op[counter][0];
+    p[1] = op[counter][1];
+    p[2] = op[counter][2];
+    m->setPoint(v, 0, p);
+    counter++;
+  }
+  m->end(vit);
+  PCU_ALWAYS_ASSERT(counter==numnp);
+
+// free op
+  for (int i = 0; i < numnp; i++)
+    free(op[i]);
+  free(op);
 }
 
 #ifdef __cplusplus