Merging develop into master

Author: Nightly Bot

README.md

@@ -53,3 +53,11 @@ For more information, start at our
 
 * If you have a usage question or have found a bug please post an [issue](https://github.com/SCOREC/core/issues).
 * Otherwise, email [Cameron Smith](https://www.scorec.rpi.edu/~cwsmith) and [email protected].
+
+### Citing PUMI
+
+If you use these tools, please cite the following paper:
+
+Daniel A. Ibanez, E. Seegyoung Seol, Cameron W. Smith, and Mark S. Shephard. 2016. PUMI: Parallel Unstructured Mesh Infrastructure. ACM Trans. Math. Softw. 42, 3, Article 17 (May 2016), 28 pages. DOI: https://doi.org/10.1145/2814935
+
+We would be happy to provide feedback on journal submissions using PUMI prior to publication.

ma/CMakeLists.txt

@@ -25,6 +25,7 @@ set(SOURCES
   maLayerCollapse.cc
   maMatch.cc
   maSolutionTransfer.cc
+  maSolutionTransferHelper.cc
   maSnap.cc
   maEdgeSwap.cc
   maShape.cc

ma/maAdapt.cc

@@ -353,6 +353,7 @@ Cavity::Cavity()
   adapter = 0;
   solutionTransfer = 0;
   shape = 0;
+  shouldTransferSizeField = false;
 }
 
 void Cavity::init(Adapt* a)
@@ -368,6 +369,8 @@ void Cavity::init(Adapt* a)
       shouldTransfer = true;
     if (shape->hasNodesOn(d))
       shouldFit = true;
+    if (adapter->sizeField->hasNodesOn(d))
+      shouldTransferSizeField = true;
   }
 }
 
@@ -409,6 +412,12 @@ void Cavity::transfer(EntityArray& oldElements)
     newEntities.retrieve(a);
     solutionTransfer->onCavity(oldElements,a);
   }
+  if (shouldTransferSizeField)
+  {
+    EntityArray a;
+    newEntities.retrieve(a);
+    adapter->sizeField->onCavity(oldElements,a);
+  }
 }
 
 void Cavity::fit(EntityArray& oldElements)

ma/maAdapt.h

@@ -137,6 +137,7 @@ class Cavity
     void fit(EntityArray& oldElements);
     bool shouldTransfer;
     bool shouldFit;
+    bool shouldTransferSizeField;
   private:
     Adapt* adapter;
     SolutionTransfer* solutionTransfer;

ma/maFaceSplit.cc

@@ -123,16 +123,22 @@ void FaceSplit::cancel()
 void FaceSplit::transfer()
 {
   Mesh* m = adapter->mesh;
+  int td;
+  td = adapter->shape->getTransferDimension();
+  for (int d = td; d <= m->getDimension(); ++d)
+    for (size_t i = 0; i < toSplit[d].getSize(); ++i)
+      adapter->shape->onRefine(toSplit[d][i], newEntities[d][i]);
+
   SolutionTransfer* st = adapter->solutionTransfer;
-  int td = st->getTransferDimension();
+  td = st->getTransferDimension();
   for (int d = td; d <= m->getDimension(); ++d)
     for (size_t i = 0; i < toSplit[d].getSize(); ++i)
       st->onRefine(toSplit[d][i], newEntities[d][i]);
 
-  td = adapter->shape->getTransferDimension();
+  td = adapter->sizeField->getTransferDimension();
   for (int d = td; d <= m->getDimension(); ++d)
     for (size_t i = 0; i < toSplit[d].getSize(); ++i)
-      adapter->shape->onRefine(toSplit[d][i], newEntities[d][i]);
+      adapter->sizeField->onRefine(toSplit[d][i], newEntities[d][i]);
 }
 
 void FaceSplit::destroyOldElements()

ma/maRefine.cc

@@ -357,6 +357,13 @@ void transferElements(Refine* r)
       for (size_t i=0; i < r->toSplit[d].getSize(); ++i)
 	st->onRefine(r->toSplit[d][i],r->newEntities[d][i]);
   }
+  {
+    SizeField* sf = a->sizeField;
+    int td = sf->getTransferDimension();
+    for (int d = td; d <= m->getDimension(); ++d)
+      for (size_t i=0; i < r->toSplit[d].getSize(); ++i)
+	sf->onRefine(r->toSplit[d][i],r->newEntities[d][i]);
+  }
 }
 
 void forgetNewEntities(Refine* r)

ma/maSize.cc

@@ -9,6 +9,7 @@
 *******************************************************************************/
 #include <PCU.h>
 #include "maSize.h"
+#include "maSolutionTransferHelper.h"
 #include "apfMatrix.h"
 #include <apfShape.h>
 #include <cstdlib>
@@ -20,6 +21,32 @@ SizeField::~SizeField()
 {
 }
 
+void SizeField::onRefine(Entity*, EntityArray&)
+{
+  PCU_ALWAYS_ASSERT_VERBOSE(0,
+      "unimplemented onRefine was called for a size-field!");
+}
+
+void SizeField::onCavity(EntityArray&, EntityArray&)
+{
+  PCU_ALWAYS_ASSERT_VERBOSE(0,
+      "unimplemented onCavity was called for a size-field!");
+}
+
+int SizeField::getTransferDimension()
+{
+  // By default there should be no size_field transfer.
+  // This is used in a loop to get the lowest dimension
+  // entities that require transfer. So something bigger
+  // than mesh dimension will ignore that loop
+  return 4;
+}
+
+bool SizeField::hasNodesOn(int)
+{
+  return false;
+}
+
 IdentitySizeField::IdentitySizeField(Mesh* m):
   mesh(m)
 {
@@ -96,14 +123,14 @@ static void orthogonalizeR(Matrix& R)
 
 static void orthogonalEigenDecompForSymmetricMatrix(Matrix const& A, Vector& v, Matrix& R)
 {
-  /* here we assume A to be real symmetric 3x3 matrix, 
+  /* here we assume A to be real symmetric 3x3 matrix,
    * we should be able to get 3 orthogonal eigen vectors
    * we also normalize the eigen vectors */
   double eigenValues[3];
   Vector eigenVectors[3];
-  
+
   apf::eigen(A, eigenVectors, eigenValues);
-  
+
   Matrix RT(eigenVectors); // eigen vectors are stored in the rows of RT
 
   RT[0] = RT[0].normalize();
@@ -132,8 +159,8 @@ static double parentMeasure[apf::Mesh::TYPES] =
 class SizeFieldIntegrator : public apf::Integrator
 {
   public:
-    SizeFieldIntegrator(SizeField* sF):
-      Integrator(2),
+    SizeFieldIntegrator(SizeField* sF, int order):
+      Integrator(order),
       measurement(0),
       sizeField(sF),
       meshElement(0),
@@ -181,7 +208,8 @@ struct MetricSizeField : public SizeField
 {
   double measure(Entity* e)
   {
-    SizeFieldIntegrator sFI(this); 
+    SizeFieldIntegrator sFI(this,
+    	std::max(mesh->getShape()->getOrder(), order)+1);
     apf::MeshElement* me = apf::createMeshElement(mesh, e);
     sFI.process(me);
     apf::destroyMeshElement(me);
@@ -201,6 +229,7 @@ struct MetricSizeField : public SizeField
     return measure(e) / parentMeasure[mesh->getType(e)];
   }
   Mesh* mesh;
+  int order; // this is the underlying sizefield order (default 1)
 };
 
 AnisotropicFunction::~AnisotropicFunction()
@@ -337,13 +366,16 @@ struct AnisoSizeField : public MetricSizeField
 {
   AnisoSizeField()
   {
+    mesh = 0;
+    order = 1;
   }
   AnisoSizeField(Mesh* m, AnisotropicFunction* f):
     bothEval(f),
     sizesEval(&bothEval),
     frameEval(&bothEval)
   {
     mesh = m;
+    order = 1;
     hField = apf::createUserField(m, "ma_sizes", apf::VECTOR,
         apf::getLagrange(1), &sizesEval);
     rField = apf::createUserField(m, "ma_frame", apf::MATRIX,
@@ -357,6 +389,7 @@ struct AnisoSizeField : public MetricSizeField
   void init(Mesh* m, apf::Field* sizes, apf::Field* frames)
   {
     mesh = m;
+    order = apf::getShape(sizes)->getOrder();
     hField = sizes;
     rField = frames;
   }
@@ -424,11 +457,14 @@ struct LogAnisoSizeField : public MetricSizeField
 {
   LogAnisoSizeField()
   {
+    mesh = 0;
+    order = 1;
   }
   LogAnisoSizeField(Mesh* m, AnisotropicFunction* f):
     logMEval(f)
   {
     mesh = m;
+    order = 1;
     logMField = apf::createUserField(m, "ma_logM", apf::MATRIX,
         apf::getLagrange(1), &logMEval);
   }
@@ -439,20 +475,34 @@ struct LogAnisoSizeField : public MetricSizeField
   void init(Mesh* m, apf::Field* sizes, apf::Field* frames)
   {
     mesh = m;
-    logMField = apf::createField(m, "ma_logM", apf::MATRIX, apf::getLagrange(1));
-    Entity* v;
-    Iterator* it = m->begin(0);
-    while ( (v = m->iterate(it)) ) {
-      Vector h;
-      Matrix f;
-      apf::getVector(sizes, v, 0, h);
-      apf::getMatrix(frames, v, 0, f);
-      Vector s(log(1/h[0]/h[0]), log(1/h[1]/h[1]), log(1/h[2]/h[2]));
-      Matrix S(s[0], 0   , 0,
-              0    , s[1], 0,
-              0    , 0   , s[2]);
-      apf::setMatrix(logMField, v, 0, f * S * transpose(f));
+    order = apf::getShape(sizes)->getOrder();
+    logMField = apf::createField(m, "ma_logM", apf::MATRIX,
+        apf::getShape(sizes));
+    int dim = m->getDimension();
+    Entity* ent;
+    Iterator* it;
+    for (int d = 0; d <= dim; d++) {
+      if (!apf::getShape(logMField)->countNodesOn(apf::Mesh::simplexTypes[d]))
+        continue;
+      it = m->begin(d);
+      while( (ent = m->iterate(it)) ){
+	int type = m->getType(ent);
+	int non = apf::getShape(logMField)->countNodesOn(type);
+	for (int i = 0; i < non; i++) {
+	  Vector h;
+	  Matrix f;
+	  apf::getVector(sizes, ent, i, h);
+	  apf::getMatrix(frames, ent, i, f);
+	  Vector s(log(1/h[0]/h[0]), log(1/h[1]/h[1]), log(1/h[2]/h[2]));
+	  Matrix S(s[0], 0   , 0,
+	      0    , s[1], 0,
+	      0    , 0   , s[2]);
+	  apf::setMatrix(logMField, ent, i, f * S * transpose(f));
+	}
+      }
+      m->end(it);
     }
+    fieldVal.allocate(apf::countComponents(logMField));
   }
   void getTransform(
       apf::MeshElement* me,
@@ -497,6 +547,34 @@ struct LogAnisoSizeField : public MetricSizeField
                           0,value,0,
                           0,0,value));
   }
+  void onRefine(
+      Entity* parent,
+      EntityArray& newEntities)
+  {
+    transfer(logMField, &(fieldVal[0]), 1, &parent, newEntities);
+  }
+  void onCavity(
+      EntityArray& oldElements,
+      EntityArray& newEntities)
+  {
+    transfer(logMField, &(fieldVal[0]),
+        oldElements.getSize(), &(oldElements[0]), newEntities);
+  }
+  int getTransferDimension()
+  {
+    int transferDimension = 4;
+    for (int d = 1; d <=3; d++)
+      if (hasNodesOn(d)) {
+        transferDimension = d;
+        break;
+      }
+    return transferDimension;
+  }
+  bool hasNodesOn(int dimension)
+  {
+    return apf::getShape(logMField)->hasNodesIn(dimension);
+  }
+  apf::NewArray<double> fieldVal;
   apf::Field* logMField;
   LogMEval logMEval;
 };

ma/maSize.h

@@ -36,6 +36,14 @@ class SizeField
         Vector const& xi,
         Matrix& t) = 0;
     virtual double getWeight(Entity* e) = 0;
+    virtual void onRefine(
+        Entity* parent,
+        EntityArray& newEntities);
+    virtual void onCavity(
+        EntityArray& oldElements,
+        EntityArray& newEntities);
+    virtual int getTransferDimension();
+    virtual bool hasNodesOn(int dimension);
 };
 
 struct IdentitySizeField : public SizeField
@@ -54,6 +62,7 @@ struct IdentitySizeField : public SizeField
           Matrix& t);
   double getWeight(Entity*);
   Mesh* mesh;
+
 };
 
 struct UniformRefiner : public IdentitySizeField