Initial implementation of support for complex fields

CMakeLists.txt

@@ -84,6 +84,13 @@ endif()
 
 set(Trilinos_PREFIX "" CACHE STRING "Trilinos installation directory")
 
+option(PUMI_CXX_COMPLEX "Use std::complex<double> instead of double_complex" OFF)
+if(PUMI_CXX_COMPLEX)
+  set(PUMI_COMPLEX PUMI_CXX_COMPLEX)
+else()
+  set(PUMI_COMPLEX PUMI_C_COMPLEX)
+endif()
+
 option(SKIP_SIMMETRIX_VERSION_CHECK "enable at your own risk; it may result in undefined behavior" OFF)
 option(ENABLE_SIMMETRIX "Build with Simmetrix support" OFF)
 message(STATUS "ENABLE_SIMMETRIX: ${ENABLE_SIMMETRIX}")

apf/CMakeLists.txt

@@ -3,6 +3,9 @@ if(DEFINED TRIBITS_PACKAGE)
   return()
 endif()
 
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/apfComplexType.h.in
+               ${CMAKE_CURRENT_SOURCE_DIR}/apfComplexType.h @ONLY)
+
 # Package sources
 set(SOURCES
   apf.cc
@@ -29,6 +32,7 @@ set(SOURCES
   apfVectorElement.cc
   apfVectorField.cc
   apfPackedField.cc
+  apfComplexField.cc
   apfNumbering.cc
   apfMixedNumbering.cc
   apfAdjReorder.cc
@@ -47,6 +51,7 @@ set(SOURCES
   apfSimplexAngleCalcs.cc
   apfFile.cc
   apfMIS.cc
+  apfZero.cc
 )
 
 # Package headers
@@ -73,6 +78,9 @@ set(HEADERS
   apfField.h
   apfFieldData.h
   apfNumberingClass.h
+  apfComplex.h
+  apfComplexType.h
+  apfElementType.h
 )
 
 # Add the apf library
@@ -92,7 +100,7 @@ target_link_libraries(apf
      lion
      can
      mth
-   )
+     )
 
 scorec_export_library(apf)
 

apf/apf.cc

@@ -6,13 +6,15 @@
  */
 
 #include "apf.h"
+#include "apfElement.h"
 #include "apfScalarField.h"
 #include "apfScalarElement.h"
 #include "apfVectorField.h"
 #include "apfVectorElement.h"
 #include "apfMatrixField.h"
 #include "apfMatrixElement.h"
 #include "apfPackedField.h"
+#include "apfComplexField.h"
 #include "apfIntegrate.h"
 #include "apfArrayData.h"
 #include "apfTagData.h"
@@ -74,7 +76,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 +128,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)
@@ -371,7 +377,7 @@ double computeCosAngle(Mesh* m, MeshEntity* pe, MeshEntity* e1, MeshEntity* e2,
 
 int countNodes(Element* e)
 {
-  return e->getShape()->countNodes();
+  return countNodes<double>(e);
 }
 
 void getScalarNodes(Element* e, NewArray<double>& values)
@@ -392,16 +398,24 @@ void getMatrixNodes(Element* e, NewArray<Matrix3x3>& values)
   element->getValues(values);
 }
 
+void getPackedNodes(Element* e, NewArray<double>& values)
+{
+  FieldBase* f = e->getFieldBase();
+  int cmps = f->countComponents();
+  ElementOf<double>* element = static_cast<ElementOf<double>*>(e);
+  element->getValues(values,cmps);
+}
+
 void getShapeValues(Element* e, Vector3 const& local,
     NewArray<double>& values)
 {
-  e->getShape()->getValues(e->getMesh(), e->getEntity(), local,values);
+  getShapeValues<double>(e,local,values);
 }
 
 void getShapeGrads(Element* e, Vector3 const& local,
     NewArray<Vector3>& grads)
 {
-  e->getGlobalGradients(local,grads);
+  getShapeGrads<double>(e,local,grads);
 }
 
 FieldShape* getShape(Field* f)
@@ -473,7 +487,7 @@ void unfreeze(Field* f)
 
 bool isFrozen(Field* f)
 {
-  return f->getData()->isFrozen();
+  return isFrozen(static_cast<FieldBase*>(f));
 }
 
 Function::~Function()

apf/apf.h

@@ -9,6 +9,7 @@
 #define APF_H
 
 #include "apfMatrix.h"
+#include "apfElementType.h"
 #include "apfNew.h"
 #include "apfDynamicArray.h"
 
@@ -31,7 +32,6 @@
 namespace apf {
 
 class Field;
-class Element;
 class Mesh;
 class MeshEntity;
 class VectorElement;
@@ -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,8 +721,21 @@ 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 * getDoubleArrayData(Field* f);
+// \note deprecated, retained for legacy compatibility until C++14 switch
+inline double* getArrayData(Field* f) { return getDoubleArrayData(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.
  */
-double* getArrayData(Field* f);
+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,46 @@ 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 * getDoubleArrayData(Field * f) { return getArrayDataT<double>(f); }
+int * getIntArrayData(Field * f) { return getArrayDataT<int>(f); }
+
+class ComplexField;
+double_complex * getComplexArrayData(ComplexField * f) { return getArrayDataT<double_complex>(reinterpret_cast<FieldBase*>(f)); }
+
+
 }

apf/apfComplex.h

@@ -0,0 +1,58 @@
+#ifndef APFCOMPLEX_H_
+#define APFCOMPLEX_H_
+
+#include "apfComplexType.h"
+#include "apfElementType.h"
+
+namespace apf
+{
+
+// forward decls for the interface
+class ComplexField;
+class Mesh;
+class FieldShape;
+class MeshEntity;
+class VectorElement;
+typedef VectorElement MeshElement;
+class Vector3;
+template <class T>
+class NewArray;
+
+ComplexField* createComplexPackedField(Mesh* m,
+                                       const char* name,
+                                       int components,
+                                       FieldShape* shape = NULL);
+
+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, Vector3 const& local, NewArray<double>& values);
+void getShapeGrads(ComplexElement* e, Vector3 const& local, NewArray<double>& grades);
+void getPackedNodes(ComplexElement* e, NewArray<double_complex>& values);
+
+}
+
+#endif