Initial implementation of support for complex fields

apf/CMakeLists.txt

@@ -48,6 +48,7 @@ set(SOURCES
   apfSimplexAngleCalcs.cc
   apfFile.cc
   apfMIS.cc
+  apfZero.cc
 )
 
 # Package headers

apf/apfComplex.h

@@ -11,4 +11,57 @@
   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

apf/apfComplexField.cc

@@ -1,21 +1,122 @@
 #include "apfComplexField.h"
+#include "apfArrayData.h"
 #include "apfElement.h"
+#include "apfTagData.h"
+#include "apfZero.h"
+#include <pcu_util.h>
 
-namespace apf {
+namespace apf
+{
+
+FieldDataOf<double_complex>* ComplexField::getData()
+{
+  return static_cast<FieldDataOf<double_complex>*>(data);
+}
 
-Element * ComplexPackedField::getElement(VectorElement * e)
+ComplexElement * ComplexPackedField::getElement(VectorElement * e)
 {
-  return new Element(this,e);
+  return new ComplexElement(this,e);
 }
 
 void ComplexPackedField::project(Field*)
 {
-  fail("ComplexPackedField::project unimplemented");
+  fail("ComplexPackedField::project is unimplemented");
+}
+
+void ComplexPackedField::axpy(double_complex, Field*)
+{
+  fail("ComplexPackedField::axpy is unimplemented");
+}
+
+ComplexField* makeComplexField(Mesh* m,
+                               const char* name,
+                               int valueType,
+                               int components,
+                               FieldShape * shape,
+                               FieldData * data)
+{
+  PCU_ALWAYS_ASSERT(!m->findComplexField(name));
+  ComplexField* f = 0;
+  if(valueType == PACKED)
+    f = new ComplexPackedField(components);
+  else
+    fail("invalid valueType in complex field construction\n");
+  f->init(name,m,shape,data);
+  m->addComplexField(f);
+  return f;
 }
 
-void ComplexPackedField::axpy(double, Field*)
+ComplexField* createComplexPackedField(Mesh* m,
+                                       const char* name,
+                                       int valueType,
+                                       int components,
+                                       FieldShape* shape)
 {
-  fail("ComplexPackedField::axpy unimplemented");
+  return makeComplexField(m,name,valueType,components,shape,
+                          new TagDataOf<double_complex>);
 }
 
+void freeze(ComplexField* f)
+{
+  if(isFrozen(f)) return;
+  f->getMesh()->hasFrozenFields = true;
+  freezeFieldData<double_complex>(f);
+}
+
+void unfreeze(ComplexField* f)
+{
+  if(isFrozen(f))
+    unfreezeFieldData<double_complex>(f);
+}
+
+bool isFrozen(ComplexField* f)
+{
+  return isFrozen(static_cast<FieldBase*>(f));
+}
+
+void zero(ComplexField* f)
+{
+  ZeroOp<double_complex> op(f);
+  op.apply(f);
+}
+
+void setComponents(ComplexField* f, MeshEntity* e, int node, double_complex const * components)
+{
+  PCU_DEBUG_ASSERT(f->getValueType() == SCALAR || f->getValueType() == PACKED);
+  f->getData()->setNodeComponents(e,node,components);
+}
+
+void getComponents(ComplexField* f, MeshEntity* e, int node, double_complex * components)
+{
+  PCU_DEBUG_ASSERT(f->getValueType() == SCALAR || f->getValueType() == PACKED);
+  f->getData()->getNodeComponents(e,node,components);
+}
+
+ComplexElement* createElement(ComplexField* f, MeshElement* e)
+{
+  return f->getElement(e);
+}
+
+ComplexElement* createElement(ComplexField* f, MeshEntity* e)
+{
+  return new ComplexElement(f,e);
+}
+
+void destroyElement(ComplexElement* e)
+{
+  delete e;
+}
+
+void getComponents(ComplexElement* e, Vector3 const& param, double_complex* components)
+{
+  e->getComponents(param,components);
+}
+
+int countNodes(ComplexElement* e)
+{
+  return countNodes(static_cast<ElementBase*>(e));
+}
+
+
+
 }

apf/apfComplexField.h

@@ -1,21 +1,45 @@
 #ifndef APFCOMPLEXFIELD_H_
 #define APFCOMPLEXFIELD_H_
 
+#include "apfElement.h"
 #include "apfField.h"
 #include "apf.h"
 
 namespace apf
 {
-class ComplexPackedField : public Field
+
+class ComplexElement : public ElementT<double_complex>
+{
+public:
+  ComplexElement(FieldBase* f, MeshEntity* e)
+    : ElementT<double_complex>(f,e)
+  { }
+  ComplexElement(FieldBase* f, VectorElement* p)
+    : ElementT<double_complex>(f,p)
+  { }
+};
+
+class ComplexField : public FieldBase
+{
+public:
+  virtual ComplexElement * getElement(VectorElement*) = 0;
+  virtual int getValueType() const = 0;
+  virtual int getScalarType() { return Mesh::COMPLEX; }
+  FieldDataOf<double_complex>* getData();
+  virtual void project(Field * frm) = 0;
+  virtual void axpy(double_complex a, Field* x) = 0;
+};
+
+class ComplexPackedField : public ComplexField
 {
 public:
-  ComplexPackedField(int c) : components(c) {}
+  ComplexPackedField(int c):components(c) {}
   virtual ~ComplexPackedField() {}
-  virtual Element * getElement(VectorElement*);
-  virtual int getValueType() const { return COMPLEX_PACKED; }
+  virtual ComplexElement * getElement(VectorElement* e);
+  virtual int getValueType() const { return PACKED; }
   virtual int countComponents() const { return components; }
-  virtual void project(Field  * frm);
-  virtual void axpy(double a, Field * x);
+  virtual void project(Field* frm);
+  virtual void axpy(double_complex a, Field* x);
 private:
   int components;
 };

apf/apfField.cc

@@ -6,8 +6,10 @@
  */
 
 #include "apfField.h"
+#include "apfComplexField.h"
 #include "apfShape.h"
 #include "apfTagData.h"
+#include "apfZero.h"
 
 namespace apf {
 
@@ -82,8 +84,8 @@ void FieldOp::apply(FieldBase* f)
     while ((e = m->iterate(it)))
     {
       // this condition makes sure that the entity has node(s)
-	  if ( s->countNodesOn(m->getType(e)) == 0)
-		continue;
+      if ( s->countNodesOn(m->getType(e)) == 0)
+        continue;
       if ( ! this->inEntity(e))
         continue;
       int n = f->countNodesOn(e);
@@ -95,34 +97,9 @@ void FieldOp::apply(FieldBase* f)
   }
 }
 
-struct ZeroOp : public FieldOp
-{
-  ZeroOp(Field* f)
-  {
-    field = f;
-    int n = f->countComponents();
-    data.allocate(n);
-    for (int i = 0; i < n; ++i)
-      data[i] = 0;
-    ent = 0;
-  }
-  bool inEntity(MeshEntity* e)
-  {
-    ent = e;
-    return true;
-  }
-  void atNode(int n)
-  {
-    setComponents(field, ent, n, &data[0]);
-  }
-  Field* field;
-  MeshEntity* ent;
-  apf::NewArray<double> data;
-};
-
 void zeroField(Field* f)
 {
-  ZeroOp op(f);
+  ZeroOp<double> op(f);
   op.apply(f);
 }
 

apf/apfField.h

@@ -71,12 +71,12 @@ class FieldOp
 };
 
 Field* makeField(
-    Mesh* m,
-    const char* name,
-    int valueType,
-    int components,
-    FieldShape* shape,
-    FieldData* data);
+  Mesh* m,
+  const char* name,
+  int valueType,
+  int components,
+  FieldShape* shape,
+  FieldData* data);
 
 bool isFrozen(FieldBase * fb);
 

apf/apfFieldData.cc

@@ -297,6 +297,8 @@ template void copyFieldData<double>(
     FieldDataOf<double>* from, FieldDataOf<double>* to);
 template void copyFieldData<long>(
     FieldDataOf<long>* from, FieldDataOf<long>* to);
+template void copyFieldData<double_complex>(
+    FieldDataOf<double_complex> * from, FieldDataOf<double_complex>* to);
 
 template <class T>
 class MultiplyOp : public FieldOp

apf/apfFieldOf.h

@@ -13,7 +13,10 @@
 
 namespace apf {
 
-template <class T>
+template <class E>
+using enable = typename std::enable_if_t<std::is_standard_layout<E>::value, bool>;
+
+template <class T, typename = void>
 class FieldOf;
 
 template <class T>
@@ -24,7 +27,7 @@ template <class T>
 void axpy(double a, FieldOf<T>* x, FieldOf<T>* y);
 
 template <class T>
-class FieldOf : public Field
+class FieldOf<T, typename std::enable_if<std::is_standard_layout<T>::value>::type > : public Field
 {
   public:
     virtual ~FieldOf() {}
@@ -49,5 +52,4 @@ class FieldOf : public Field
 };
 
 }
-
 #endif