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
@@ -73,6 +74,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"
@@ -160,13 +161,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 +178,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 +195,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 +212,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);
 }
 
@@ -473,7 +484,7 @@ void unfreeze(Field* f)
 
 bool isFrozen(Field* f)
 {
-  return f->getData()->isFrozen();
+  return isFrozen(static_cast<FieldBase*>(f));
 }
 
 Function::~Function()

apf/apf.h

@@ -11,6 +11,7 @@
 #include "apfMatrix.h"
 #include "apfNew.h"
 #include "apfDynamicArray.h"
+#include "apfComplex.h"
 
 #include <vector>
 #include <limits>
@@ -717,9 +718,29 @@ 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 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(Field * f);
+
 /** \brief Initialize all nodal values with all-zero components */
 void zeroField(Field* f);
 

apf/apfArrayData.cc

@@ -1,4 +1,5 @@
 #include "apfArrayData.h"
+#include "apfComplex.h"
 #include "apfNumbering.h"
 #include "apfTagData.h"
 
@@ -125,19 +126,29 @@ 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)
+{
+  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);
+
+
 }

apf/apfComplex.h

@@ -0,0 +1,14 @@
+#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
+
+#endif

apf/apfComplexField.cc

@@ -0,0 +1,21 @@
+#include "apfComplexField.h"
+#include "apfElement.h"
+
+namespace apf {
+
+Element * ComplexPackedField::getElement(VectorElement * e)
+{
+  return new Element(this,e);
+}
+
+void ComplexPackedField::project(Field*)
+{
+  fail("ComplexPackedField::project unimplemented");
+}
+
+void ComplexPackedField::axpy(double, Field*)
+{
+  fail("ComplexPackedField::axpy unimplemented");
+}
+
+}

apf/apfComplexField.h

@@ -0,0 +1,25 @@
+#ifndef APFCOMPLEXFIELD_H_
+#define APFCOMPLEXFIELD_H_
+
+#include "apfField.h"
+#include "apf.h"
+
+namespace apf
+{
+class ComplexPackedField : public Field
+{
+public:
+  ComplexPackedField(int c) : components(c) {}
+  virtual ~ComplexPackedField() {}
+  virtual Element * getElement(VectorElement*);
+  virtual int getValueType() const { return COMPLEX_PACKED; }
+  virtual int countComponents() const { return components; }
+  virtual void project(Field  * frm);
+  virtual void axpy(double a, Field * x);
+private:
+  int components;
+};
+
+}
+
+#endif

apf/apfField.cc

@@ -126,4 +126,9 @@ void zeroField(Field* f)
   op.apply(f);
 }
 
+bool isFrozen(FieldBase * fb)
+{
+  return fb->getData()->isFrozen();
+}
+
 } //namespace apf

apf/apfField.h

@@ -78,6 +78,8 @@ Field* makeField(
     FieldShape* shape,
     FieldData* data);
 
+bool isFrozen(FieldBase * fb);
+
 } //namespace apf
 
 #endif

apf/apfFieldData.cc

@@ -1,4 +1,5 @@
 #include <PCU.h>
+#include "apfComplex.h"
 #include "apfFieldData.h"
 #include "apfShape.h"
 #include <pcu_util.h>
@@ -72,6 +73,7 @@ void synchronizeFieldData(FieldDataOf<T>* data, Sharing* shr, bool delete_shr)
 }
 
 /* instantiate here */
+template void synchronizeFieldData<double_complex>(FieldDataOf<double_complex>*, Sharing*, bool);
 template void synchronizeFieldData<int>(FieldDataOf<int>*, Sharing*, bool);
 template void synchronizeFieldData<double>(FieldDataOf<double>*, Sharing*, bool);
 template void synchronizeFieldData<long>(FieldDataOf<long>*, Sharing*, bool);
@@ -250,6 +252,7 @@ int FieldDataOf<T>::getElementData(MeshEntity* entity, NewArray<T>& data)
   return n;
 }
 
+template class FieldDataOf<double_complex>;
 template class FieldDataOf<double>;
 template class FieldDataOf<int>;
 template class FieldDataOf<long>;

apf/apfMesh.h

@@ -14,6 +14,7 @@
 #include <vector>
 #include <map>
 #include <set>
+#include "apfComplex.h"
 #include "apfVector.h"
 #include "apfDynamicArray.h"
 
@@ -225,6 +226,8 @@ class Mesh
     virtual void getResidence(MeshEntity* e, Parts& residence) = 0;
     /** \brief Creates a double array tag over the mesh given a name and size */
     virtual MeshTag* createDoubleTag(const char* name, int size) = 0;
+    /** \brief creates a double_complex array tag over the mesh given a name and size */
+    virtual MeshTag* createComplexTag(const char* name, int size) = 0;
     /** \brief Creates an int array tag over the mesh given a name and size */
     virtual MeshTag* createIntTag(const char* name, int size) = 0;
     /** \brief Creates a long array tag over the mesh given a name and size */
@@ -239,6 +242,10 @@ class Mesh
     virtual void getDoubleTag(MeshEntity* e, MeshTag* tag, double* data) = 0;
     /** \brief set double array tag data */
     virtual void setDoubleTag(MeshEntity* e, MeshTag* tag, double const* data) = 0;
+    /** \brief get a complex array tag */
+    virtual void getComplexTag(MeshEntity* e, MeshTag* tag, double_complex* data) = 0;
+    /** \brief set a complex array tag */
+    virtual void setComplexTag(MeshEntity* e, MeshTag* tag, double_complex const * data) = 0;
     /** \brief get int array tag data */
     virtual void getIntTag(MeshEntity* e, MeshTag* tag, int* data) = 0;
     /** \brief set int array tag data */
@@ -263,7 +270,9 @@ class Mesh
       /** \brief signed 32-bit integer */
       INT,
       /** \brief signed 64-bit integer */
-      LONG };
+      LONG,
+      /** \brief either std::complex<double> or double_complex (C) depending on configuration */
+      COMPLEX};
     /** \brief get the data type of a tag
         \return a value in apf::Mesh::TagType */
     virtual int getTagType(MeshTag* t) = 0;

apf/apfTagData.h

@@ -57,6 +57,24 @@ class TagHelper<double> : public TagMaker
     }
 };
 
+template <>
+class TagHelper<double_complex> : public TagMaker
+{
+public:
+  MeshTag * make(Mesh * m, const char * n, int s)
+  {
+    return m->createComplexTag(n,s);
+  }
+  void get(Mesh * m, MeshEntity * e, MeshTag * t, double_complex * data)
+  {
+    m->getComplexTag(e,t,data);
+  }
+  void set(Mesh * m, MeshEntity * e, MeshTag * t, double_complex const* data)
+  {
+    return m->setComplexTag(e,t,data);
+  }
+};
+
 template <>
 class TagHelper<long> : public TagMaker
 {

mds/apfMDS.cc

@@ -407,6 +407,13 @@ class MeshMDS : public Mesh2
           sizeof(double)*size, Mesh::DOUBLE);
       return reinterpret_cast<MeshTag*>(tag);
     }
+    MeshTag* createComplexTag(const char* name, int size)
+    {
+      mds_tag* tag;
+      PCU_ALWAYS_ASSERT(!mds_find_tag(&mesh->tags,name));
+      tag = mds_create_tag(&(mesh->tags),name,sizeof(double_complex)*size, Mesh::COMPLEX);
+      return reinterpret_cast<MeshTag*>(tag);
+    }
     MeshTag* createIntTag(const char* name, int size)
     {
       mds_tag* tag;
@@ -474,6 +481,14 @@ class MeshMDS : public Mesh2
     {
       setTag(e,tag,data);
     }
+    void getComplexTag(MeshEntity* e, MeshTag * tag, double_complex * data)
+    {
+      getTag(e,tag,data);
+    }
+    void setComplexTag(MeshEntity* e, MeshTag * tag, double_complex const* data)
+    {
+      setTag(e,tag,data);
+    }
     void getIntTag(MeshEntity* e, MeshTag* tag, int* data)
     {
       getTag(e,tag,data);

pumi/pumi_ghost.cc

@@ -400,7 +400,7 @@ void pumi_ghost_create(pMesh m, Ghosting* plan)
   std::vector<apf::Field*> frozen_fields;
   for (int i=0; i<m->countFields(); ++i)
   {
-    pField f = m->getField(i);
+    apf::Field * f = m->getField(i);
     if (isFrozen(f))
     {
       frozen_fields.push_back(f); // turn field data from tag to array