Initial implementation of support for complex fields

apf/CMakeLists.txt

@@ -29,6 +29,7 @@ set(SOURCES
   apfVectorElement.cc
   apfVectorField.cc
   apfPackedField.cc
+  apfComplexField.cc
   apfNumbering.cc
   apfMixedNumbering.cc
   apfAdjReorder.cc
@@ -47,6 +48,7 @@ set(SOURCES
   apfSimplexAngleCalcs.cc
   apfFile.cc
   apfMIS.cc
+  apfZero.cc
 )
 
 # Package headers
@@ -73,6 +75,7 @@ set(HEADERS
   apfField.h
   apfFieldData.h
   apfNumberingClass.h
+  apfComplex.h
 )
 
 # Add the apf library

apf/apf.cc

@@ -13,6 +13,7 @@
 #include "apfMatrixField.h"
 #include "apfMatrixElement.h"
 #include "apfPackedField.h"
+#include "apfComplexField.h"
 #include "apfIntegrate.h"
 #include "apfArrayData.h"
 #include "apfTagData.h"
@@ -74,7 +75,7 @@ Field* makeField(
   else if (valueType == PACKED)
     f = new PackedField(components);
   else
-    fail("invalid valueType in field construction\n");
+    fail("invalid valueType in double field construction\n");
   f->init(name,m,shape,data);
   m->addField(f);
   return f;
@@ -126,8 +127,12 @@ Field* createPackedField(Mesh* m, const char* name, int components,
 
 Field* cloneField(Field* f, Mesh* onto)
 {
-  return makeField(onto, f->getName(), f->getValueType(), f->countComponents(),
-      f->getShape(), f->getData()->clone());
+  return makeField(onto,
+                   f->getName(),
+                   f->getValueType(),
+                   f->countComponents(),
+                   f->getShape(),
+                   f->getData()->clone());
 }
 
 Mesh* getMesh(Field* f)
@@ -160,13 +165,15 @@ void destroyField(Field* f)
 
 void setScalar(Field* f, MeshEntity* e, int node, double value)
 {
+  PCU_DEBUG_ASSERT(f->getValueType() == SCALAR);
   ScalarField* field = static_cast<ScalarField*>(f);
   double tmp[1] = {value};
   field->setNodeValue(e,node,tmp);
 }
 
 double getScalar(Field* f, MeshEntity* e, int node)
 {
+  PCU_DEBUG_ASSERT(f->getValueType() == SCALAR);
   ScalarField* field = static_cast<ScalarField*>(f);
   double value[1];
   field->getNodeValue(e,node,value);
@@ -175,13 +182,15 @@ double getScalar(Field* f, MeshEntity* e, int node)
 
 void setVector(Field* f, MeshEntity* e, int node, Vector3 const& value)
 {
+  PCU_DEBUG_ASSERT(f->getValueType() == VECTOR);
   VectorField* field = static_cast<VectorField*>(f);
   Vector3 tmp[1] = {value};
   field->setNodeValue(e,node,tmp);
 }
 
 void getVector(Field* f, MeshEntity* e, int node, Vector3& value)
 {
+  PCU_DEBUG_ASSERT(f->getValueType() == VECTOR);
   VectorField* field = static_cast<VectorField*>(f);
   Vector3 tmp[1];
   field->getNodeValue(e,node,tmp);
@@ -190,13 +199,15 @@ void getVector(Field* f, MeshEntity* e, int node, Vector3& value)
 
 void setMatrix(Field* f, MeshEntity* e, int node, Matrix3x3 const& value)
 {
+  PCU_DEBUG_ASSERT(f->getValueType() == MATRIX);
   MatrixField* field = static_cast<MatrixField*>(f);
   Matrix3x3 tmp[1] = {value};
   field->setNodeValue(e,node,tmp);
 }
 
 void getMatrix(Field* f, MeshEntity* e, int node, Matrix3x3& value)
 {
+  PCU_DEBUG_ASSERT(f->getValueType() == MATRIX);
   MatrixField* field = static_cast<MatrixField*>(f);
   Matrix3x3 tmp[1];
   field->getNodeValue(e,node,tmp);
@@ -205,11 +216,15 @@ void getMatrix(Field* f, MeshEntity* e, int node, Matrix3x3& value)
 
 void setComponents(Field* f, MeshEntity* e, int node, double const* components)
 {
+  PCU_DEBUG_ASSERT(f->getScalarType() == Mesh::DOUBLE &&
+                   (f->getValueType() == SCALAR || f->getValueType() == PACKED));
   f->getData()->setNodeComponents(e,node,components);
 }
 
 void getComponents(Field* f, MeshEntity* e, int node, double* components)
 {
+  PCU_DEBUG_ASSERT(f->getScalarType() == Mesh::DOUBLE &&
+                   (f->getValueType() == SCALAR || f->getValueType() == PACKED));
   f->getData()->getNodeComponents(e,node,components);
 }
 
@@ -240,49 +255,57 @@ MeshEntity* getMeshEntity(Element* e)
 
 double getScalar(Element* e, Vector3 const& param)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<ScalarElement*>(e) != NULL);
   ScalarElement* element = static_cast<ScalarElement*>(e);
   return element->getValue(param);
 }
 
 void getGrad(Element* e, Vector3 const& param, Vector3& g)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<ScalarElement*>(e) != NULL);
   ScalarElement* element = static_cast<ScalarElement*>(e);
   element->grad(param,g);
 }
 
 void getVector(Element* e, Vector3 const& param, Vector3& value)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<VectorElement*>(e) != NULL);
   VectorElement* element = static_cast<VectorElement*>(e);
   value = element->getValue(param);
 }
 
 double getDiv(Element* e, Vector3 const& param)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<VectorElement*>(e) != NULL);
   VectorElement* element = static_cast<VectorElement*>(e);
   return element->div(param);
 }
 
 void getCurl(Element* e, Vector3 const& param, Vector3& c)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<VectorElement*>(e) != NULL);
   VectorElement* element = static_cast<VectorElement*>(e);
   return element->curl(param,c);
 }
 
 void getVectorGrad(Element* e, Vector3 const& param, Matrix3x3& deriv)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<VectorElement*>(e) != NULL);
   VectorElement* element = static_cast<VectorElement*>(e);
   return element->grad(param,deriv);
 }
 
 void getMatrix(Element* e, Vector3 const& param, Matrix3x3& value)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<MatrixElement*>(e) != NULL);
   MatrixElement* element = static_cast<MatrixElement*>(e);
   value = element->getValue(param);
 }
 
 
 void getMatrixGrad(Element* e, Vector3 const& param, Vector<27>& deriv)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<MatrixElement*>(e) != NULL);
   MatrixElement* element = static_cast<MatrixElement*>(e);
   return element->grad(param,deriv);
 }
@@ -371,27 +394,39 @@ double computeCosAngle(Mesh* m, MeshEntity* pe, MeshEntity* e1, MeshEntity* e2,
 
 int countNodes(Element* e)
 {
-  return e->getShape()->countNodes();
+  return countNodes(static_cast<ElementBase*>(e));
 }
 
 void getScalarNodes(Element* e, NewArray<double>& values)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<ElementOf<double>*>(e) != NULL);
   ElementOf<double>* element = static_cast<ElementOf<double>*>(e);
   element->getValues(values);
 }
 
 void getVectorNodes(Element* e, NewArray<Vector3>& values)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<ElementOf<Vector3>*>(e) != NULL);
   ElementOf<Vector3>* element = static_cast<ElementOf<Vector3>*>(e);
   element->getValues(values);
 }
 
 void getMatrixNodes(Element* e, NewArray<Matrix3x3>& values)
 {
+  PCU_DEBUG_ASSERT(dynamic_cast<ElementOf<Matrix3x3>*>(e) != NULL);
   ElementOf<Matrix3x3>* element = static_cast<ElementOf<Matrix3x3>*>(e);
   element->getValues(values);
 }
 
+void getPackedNodes(Element* e, NewArray<double>& values)
+{
+  FieldBase* f = e->getFieldBase();
+  int cmps = f->countComponents();
+  PCU_DEBUG_ASSERT(dynamic_cast<ElementOf<double>*>(e) != NULL);
+  ElementOf<double>* element = static_cast<ElementOf<double>*>(e);
+  element->getValues(values,cmps);
+}
+
 void getShapeValues(Element* e, Vector3 const& local,
     NewArray<double>& values)
 {
@@ -473,7 +508,7 @@ void unfreeze(Field* f)
 
 bool isFrozen(Field* f)
 {
-  return f->getData()->isFrozen();
+  return isFrozen(static_cast<FieldBase*>(f));
 }
 
 Function::~Function()

apf/apf.h

@@ -43,7 +43,7 @@ template <class T> class ReductionOp;
 template <class T> class ReductionSum;
 
 /** \brief Base class for applying operations to make a Field consistent
-  * in parallel 
+  * in parallel
   * \details This function gets applied pairwise to the Field values
   * from every partition, resulting in a single unique value.  No guarantees
   * are made about the order in which this function is applied to the
@@ -569,6 +569,10 @@ void getVectorNodes(Element* e, NewArray<Vector3>& values);
   */
 void getMatrixNodes(Element* e, NewArray<Matrix3x3>& values);
 
+/** \brief Returns the element nodal values for a packed field
+  */
+void getPackedNodes(Element* e, NewArray<double>& values);
+
 /** \brief Returns the shape function values at a point
   */
 void getShapeValues(Element* e, Vector3 const& local,
@@ -717,9 +721,20 @@ bool isFrozen(Field* f);
   \details This function is only defined for fields
   which are using array storage, for which apf::isFrozen
   returns true.
+  \note If the underlying field data type is NOT double,
+  this will cause an assert fail in all compile modes.
  */
 double* getArrayData(Field* f);
 
+/** \brief Return the contiguous array storing this field.
+  \details This function is only defined for fields
+  which are using array storage, for which apf::isFrozen
+  returns true.
+  \note If the underlying field data type is NOT int,
+  this will cause an assert fail in all compile modes.
+ */
+int* getIntArrayData(Field* f);
+
 /** \brief Initialize all nodal values with all-zero components */
 void zeroField(Field* f);
 

apf/apfArrayData.cc

@@ -1,6 +1,9 @@
 #include "apfArrayData.h"
+#include "apfComplex.h"
 #include "apfNumbering.h"
 #include "apfTagData.h"
+#include <pcu_util.h>
+#include <type_traits>
 
 namespace apf {
 
@@ -125,19 +128,43 @@ void unfreezeFieldData(FieldBase* field) {
 }
 
 /* instantiate here */
+template void freezeFieldData<double_complex>(FieldBase* field);
 template void freezeFieldData<int>(FieldBase* field);
 template void freezeFieldData<double>(FieldBase* field);
+template void unfreezeFieldData<double_complex>(FieldBase * field);
 template void unfreezeFieldData<int>(FieldBase* field);
 template void unfreezeFieldData<double>(FieldBase* field);
 
-double* getArrayData(Field* f) {
-  if (!isFrozen(f)) {
+template <typename T>
+T* getArrayDataT(FieldBase* f)
+{
+  int scalar = f->getScalarType();
+  // having to assert this is terrible and if we add more field types
+  // unsustainable and bad practice, but the current other option is
+  // changing the API and being more explicit about type storage
+  // since Field assumes it has Scalars of type double
+  PCU_ALWAYS_ASSERT(
+    (scalar == Mesh::DOUBLE && std::is_same<T,double>::value) ||
+    (scalar == Mesh::INT && std::is_same<T,int>::value) ||
+    (scalar == Mesh::LONG && std::is_same<T,long>::value) ||
+    (scalar == Mesh::COMPLEX && std::is_same<T,double_complex>::value)
+    );
+  if(!isFrozen(f))
     return 0;
-  } else {
-    FieldDataOf<double>* p = f->getData();
-    ArrayDataOf<double>* a = static_cast<ArrayDataOf<double>* > (p);
+  else
+  {
+    FieldDataOf<T>* p = reinterpret_cast<FieldDataOf<T>*>(f->getData());
+    ArrayDataOf<T>* a = static_cast<ArrayDataOf<T>*>(p);
     return a->getDataArray();
   }
 }
 
+template double_complex* getArrayDataT(FieldBase* field);
+template int* getArrayDataT(FieldBase* field);
+template double* getArrayDataT(FieldBase* field);
+
+double * getArrayData(Field * f) { return getArrayDataT<double>(f); }
+int * getIntArrayData(Field * f) { return getArrayDataT<int>(f); }
+double_complex * getComplexArrayData(Field * f) { return getArrayDataT<double_complex>(f); }
+
 }

apf/apfComplex.h

@@ -0,0 +1,67 @@
+#ifndef APFCOMPLEX_H_
+#define APFCOMPLEX_H_
+
+#ifdef C_COMPLEX
+  #include <complex.h>
+#endif
+
+#define CXX_COMPLEX 1
+#ifdef CXX_COMPLEX
+  #include <complex>
+  using double_complex = std::complex<double>;
+#endif
+
+namespace apf
+{
+
+// forward decls for the interface
+class ComplexElement;
+class ComplexField;
+class Mesh;
+class FieldShape;
+class MeshEntity;
+class VectorElement;
+typedef VectorElement MeshElement;
+class Vector3;
+template <class T>
+class NewArray;
+
+ComplexField* createComplexField(Mesh* m,
+                                 const char* name,
+                                 int valueType,
+                                 int components,
+                                 FieldShape* shape);
+
+void freeze(ComplexField* f);
+void unfreeze(ComplexField* f);
+bool isFrozen(ComplexField* f);
+
+/** \brief Return the contiguous array storing this field.
+  \details This function is only defined for fields
+  which are using array storage, for which apf::isFrozen
+  returns true.
+  \note If the underlying field data type is NOT double_complex,
+  this will cause an assert fail in all compile modes.
+ */
+double_complex* getComplexArrayData(ComplexField * f);
+void zeroField(ComplexField* f);
+
+void setComponents(ComplexField* f, MeshEntity* e, int node, double_complex const * components);
+void getComponents(ComplexField* f, MeshEntity* e, int node, double_complex * components);
+
+ComplexElement* createElement(ComplexField* f, MeshElement* e);
+ComplexElement* createElement(ComplexField* f, MeshEntity* e);
+void destroyElement(ComplexElement* e);
+
+MeshElement* getMeshElement(ComplexElement* e);
+MeshEntity* getMeshEntity(ComplexElement* e);
+
+void getComponents(ComplexElement* e, Vector3 const& param, double_complex* components);
+int countNodes(ComplexElement* e);
+void getShapeValues(ComplexElement* e);
+void getShapeGrads(ComplexElement* e);
+void getPackedNodes(ComplexElement* e, NewArray<double_complex>& values);
+
+}
+
+#endif