diff --git a/CMakeLists.txt b/CMakeLists.txt
index ae9496d6c..f9ea9d1d6 100644
--- a/CMakeLists.txt
+++ b/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}")
diff --git a/apf/CMakeLists.txt b/apf/CMakeLists.txt
index 881a9fe65..fa2296936 100644
--- a/apf/CMakeLists.txt
+++ b/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)
 
diff --git a/apf/apf.cc b/apf/apf.cc
index c7b60a1a2..94eadcf88 100644
--- a/apf/apf.cc
+++ b/apf/apf.cc
@@ -6,6 +6,7 @@
  */
 
 #include "apf.h"
+#include "apfElement.h"
 #include "apfScalarField.h"
 #include "apfScalarElement.h"
 #include "apfVectorField.h"
@@ -13,6 +14,7 @@
 #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()
diff --git a/apf/apf.h b/apf/apf.h
index 1ea32d2df..977885b02 100644
--- a/apf/apf.h
+++ b/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);
diff --git a/apf/apfArrayData.cc b/apf/apfArrayData.cc
index 3d76f7c0a..0cf365542 100644
--- a/apf/apfArrayData.cc
+++ b/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)); }
+
+
 }
diff --git a/apf/apfComplex.h b/apf/apfComplex.h
new file mode 100644
index 000000000..78301b253
--- /dev/null
+++ b/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
diff --git a/apf/apfComplexField.cc b/apf/apfComplexField.cc
new file mode 100644
index 000000000..457ee5d7d
--- /dev/null
+++ b/apf/apfComplexField.cc
@@ -0,0 +1,145 @@
+#include "apfComplexField.h"
+#include "apfArrayData.h"
+#include "apfElementOf.h"
+#include "apfTagData.h"
+#include "apfZero.h"
+#include <pcu_util.h>
+
+namespace apf
+{
+
+FieldDataOf<double_complex>* ComplexField::getData()
+{
+  return static_cast<FieldDataOf<double_complex>*>(data);
+}
+
+ComplexElement * ComplexPackedField::getElement(VectorElement * e)
+{
+  return new ComplexElement(this,e);
+}
+
+void ComplexPackedField::project(Field*)
+{
+  fail("ComplexPackedField::project is unimplemented");
+}
+
+void ComplexPackedField::axpy(double_complex, Field*)
+{
+  fail("ComplexPackedField::axpy is unimplemented");
+}
+
+ComplexField* makeComplexGeneralField(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;
+}
+
+ComplexField* createComplexPackedField(Mesh* m,
+                                       const char* name,
+                                       int components,
+                                       FieldShape* shape)
+{
+  if(shape == NULL)
+    shape = m->getShape();
+  return makeComplexGeneralField(m,name,PACKED,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 zeroField(ComplexField* f)
+{
+  ZeroOp<double_complex> op(f);
+  op.apply(f);
+}
+
+void setComponents(ComplexField* f, MeshEntity* e, int node, double_complex const * components)
+{
+  f->getData()->setNodeComponents(e,node,components);
+}
+
+void getComponents(ComplexField* f, MeshEntity* e, int node, double_complex * components)
+{
+  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;
+}
+
+MeshElement* getMeshElement(ComplexElement* e)
+{
+  return e->getParent();
+}
+
+MeshEntity* getMeshEntity(ComplexElement* e)
+{
+  return e->getEntity();
+}
+
+void getComponents(ComplexElement* e, Vector3 const& param, double_complex* components)
+{
+  e->getComponents(param,components);
+}
+
+int countNodes(ComplexElement* e)
+{
+  return countNodes<double_complex>(e);
+}
+
+void getShapeValues(ComplexElement* e, Vector3 const& local, NewArray<double>& values)
+{
+  getShapeValues<double_complex>(e,local,values);
+}
+void getShapeGrads(ComplexElement* e, Vector3 const& local, NewArray<Vector3>& grads)
+{
+  getShapeGrads<double_complex>(e,local,grads);
+}
+void getPackedNodes(ComplexElement* e, NewArray<double_complex>& values)
+{
+  FieldBase * f = e->getFieldBase();
+  int cmps = f->countComponents();
+  ComplexElementOf<double_complex>* element = static_cast<ComplexElementOf<double_complex>*>(e);
+  element->getValues(values,cmps);
+}
+
+}
diff --git a/apf/apfComplexField.h b/apf/apfComplexField.h
new file mode 100644
index 000000000..2c22ca8c4
--- /dev/null
+++ b/apf/apfComplexField.h
@@ -0,0 +1,102 @@
+#ifndef APFCOMPLEXFIELD_H_
+#define APFCOMPLEXFIELD_H_
+
+#include "apfElement.h"
+#include "apfField.h"
+#include "apf.h"
+
+namespace apf
+{
+
+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;
+};
+
+template <class T, typename = void>
+class ComplexFieldOf;
+
+template <class T>
+void project(ComplexFieldOf<T>* to, ComplexFieldOf<T>* from) {}
+
+template <class T>
+void axpy(double_complex a, ComplexFieldOf<T>* x, ComplexFieldOf<T>* y) {}
+
+template <class T>
+class ComplexFieldOf<T,typename std::enable_if<std::is_standard_layout<T>::value>::type > : public ComplexField
+{
+public:
+  virtual ~ComplexFieldOf() {}
+  void setNodeValue(MeshEntity* e, int node, T const value)
+  {
+    getData()->setNodeComponents(
+      e,node,reinterpret_cast<double_complex const*>(value));
+  }
+  void getNodeValue(MeshEntity* e, int node, T* value)
+  {
+    getData()->getNodeComponents(
+      e,node,reinterpret_cast<double_complex*>(value));
+  }
+  void project(ComplexField* from)
+  {
+    apf::project<T>(this,static_cast<ComplexFieldOf<T>*>(from));
+  }
+  void axpy(double_complex a, ComplexField* x)
+  {
+    apf::axpy<T>(this,a,reinterpret_cast<ComplexFieldOf<T>*>(x));
+  }
+};
+
+template <class T>
+class ComplexElementOf : public ComplexElement
+{
+public:
+  ComplexElementOf(ComplexFieldOf<T>* f, MeshEntity* e)
+    : Element(f,e)
+  { }
+  ComplexElementOf(ComplexFieldOf<T>* f, VectorElement* p)
+    : Element(f,p)
+  { }
+  virtual ~ComplexElementOf() { }
+  T* getNodeValues()
+  {
+    return reinterpret_cast<T*>(&(this->nodeData[0]));
+  }
+  T getValue(Vector3 const& local)
+  {
+    T value[1];
+    getComponents(local, reinterpret_cast<double_complex*>(value));
+    return value[0];
+  }
+  void getValues(NewArray<T>& values, int nc = 1)
+  {
+    values.allocate(nen * nc);
+    T* nodeValues = getNodeValues();
+    for(int i=0; i < nen * nc; ++i)
+      values[i] = nodeValues[i];
+  }
+};
+
+class ComplexPackedField : public ComplexField
+{
+public:
+  ComplexPackedField(int c):components(c) {}
+  virtual ~ComplexPackedField() {}
+  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_complex a, Field* x);
+private:
+  int components;
+};
+
+}
+
+#endif
diff --git a/apf/apfComplexType.h.in b/apf/apfComplexType.h.in
new file mode 100644
index 000000000..6ea989de1
--- /dev/null
+++ b/apf/apfComplexType.h.in
@@ -0,0 +1,10 @@
+#define @PUMI_COMPLEX@
+
+#ifndef PUMI_CXX_COMPLEX
+#include <complex.h>
+typedef double _Complex double_complex;
+#else
+  #include <complex>
+  using double_complex = std::complex<double>;
+#endif
+
diff --git a/apf/apfElement.cc b/apf/apfElement.cc
index 82c89c602..70e1e8c67 100644
--- a/apf/apfElement.cc
+++ b/apf/apfElement.cc
@@ -8,36 +8,39 @@
 #include "apfElement.h"
 #include "apfShape.h"
 #include "apfMesh.h"
+#include "apfComplexType.h"
 #include "apfVectorElement.h"
 
 namespace apf {
 
-void Element::init(Field* f, MeshEntity* e, VectorElement* p)
-{
-  field = f;
-  mesh = f->getMesh();
-  entity = e;
-  shape = f->getShape()->getEntityShape(mesh->getType(e));
-  parent = p;
-  nen = shape->countNodes();
-  nc = f->countComponents();
-  getNodeData();
-}
-
-Element::Element(Field* f, MeshEntity* e)
+template <class T>
+ElementBase<T>::ElementBase(FieldBase* f, MeshEntity* e)
 {
   init(f,e,0);
 }
 
-Element::Element(Field* f, VectorElement* p)
+template <class T>
+ElementBase<T>::ElementBase(FieldBase* f, VectorElement* p)
 {
   init(f,p->getEntity(),p);
 }
 
-Element::~Element()
+template <class T>
+void ElementBase<T>::getGlobalGradients(Vector3 const& local, NewArray<Vector3>& globalGradients)
 {
+  Matrix3x3 J;
+  parent->getJacobian(local,J);
+  Matrix3x3 jinv = getJacobianInverse(J, getDimension());
+  NewArray<Vector3> localGradients;
+  shape->getLocalGradients(mesh, entity, local,localGradients);
+  globalGradients.allocate(nen);
+  for (int i=0; i < nen; ++i)
+    globalGradients[i] = jinv * localGradients[i];
 }
 
+template class ElementBase<double>;
+template class ElementBase<double_complex>;
+
 Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim)
 {
   switch (dim) {
@@ -72,33 +75,4 @@ Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim)
   }
 }
 
-void Element::getGlobalGradients(Vector3 const& local,
-                                 NewArray<Vector3>& globalGradients)
-{
-  Matrix3x3 J;
-  parent->getJacobian(local,J);
-  Matrix3x3 jinv = getJacobianInverse(J, getDimension());
-  NewArray<Vector3> localGradients;
-  shape->getLocalGradients(mesh, entity, local,localGradients);
-  globalGradients.allocate(nen);
-  for (int i=0; i < nen; ++i)
-    globalGradients[i] = jinv * localGradients[i];
-}
-
-void Element::getComponents(Vector3 const& xi, double* c)
-{
-  NewArray<double> shapeValues;
-  shape->getValues(mesh, entity, xi, shapeValues);
-  for (int ci = 0; ci < nc; ++ci)
-    c[ci] = 0;
-  for (int ni = 0; ni < nen; ++ni)
-    for (int ci = 0; ci < nc; ++ci)
-      c[ci] += nodeData[ni * nc + ci] * shapeValues[ni];
-}
-
-void Element::getNodeData()
-{
-  field->getData()->getElementData(entity,nodeData);
-}
-
 }//namespace apf
diff --git a/apf/apfElement.h b/apf/apfElement.h
index f100d42e8..1dc17d424 100644
--- a/apf/apfElement.h
+++ b/apf/apfElement.h
@@ -10,43 +10,87 @@
 
 #include "apfMesh.h"
 #include "apfField.h"
+#include "apfFieldData.h"
 #include "apfShape.h"
 
 namespace apf {
 
+Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim);
+
 class EntityShape;
 class VectorElement;
 
-class Element
+template <class T>
+class ElementBase
 {
-  public:
-    Element(Field* f, MeshEntity* e);
-    Element(Field* f, VectorElement* p);
-    virtual ~Element();
-    void getGlobalGradients(Vector3 const& local,
-                            NewArray<Vector3>& globalGradients);
-    int getType() {return mesh->getType(entity);}
-    int getDimension() {return Mesh::typeDimension[getType()];}
-    int getOrder() {return field->getShape()->getOrder();}
-    VectorElement* getParent() {return parent;}
-    MeshEntity* getEntity() {return entity;}
-    Mesh* getMesh() {return mesh;}
-    EntityShape* getShape() {return shape;}
-    void getComponents(Vector3 const& xi, double* c);
-  protected:
-    void init(Field* f, MeshEntity* e, VectorElement* p);
-    void getNodeData();
-    Field* field;
-    Mesh* mesh;
-    MeshEntity* entity;
-    EntityShape* shape;
-    VectorElement* parent;
-    int nen;
-    int nc;
-    NewArray<double> nodeData;
+public:
+  ElementBase(FieldBase* f, MeshEntity* e);
+  ElementBase(FieldBase* f, VectorElement* p);
+  virtual ~ElementBase() {}
+  void getGlobalGradients(Vector3 const& local,
+                          NewArray<Vector3>& globalGradients);
+  int getType() { return mesh->getType(entity); }
+  int getDimension() { return Mesh::typeDimension[getType()]; }
+  int getOrder() { return field->getShape()->getOrder(); }
+  VectorElement* getParent() { return parent; }
+  MeshEntity* getEntity() { return entity; }
+  Mesh* getMesh() { return mesh; }
+  EntityShape * getShape() { return shape; }
+  FieldBase* getFieldBase() { return field; }
+  void getComponents(Vector3 const& xi, T * c)
+  {
+    NewArray<double> shapeValues;
+    shape->getValues(mesh, entity, xi, shapeValues);
+    for (int ci = 0; ci < nc; ++ci)
+      c[ci] = 0;
+    for (int ni = 0; ni < nen; ++ni)
+      for (int ci = 0; ci < nc; ++ci)
+        c[ci] += nodeData[ni * nc + ci] * shapeValues[ni];
+  }
+protected:
+  void init(FieldBase* f, MeshEntity* e, VectorElement* p)
+  {
+    field = f;
+    mesh = f->getMesh();
+    entity = e;
+    shape = f->getShape()->getEntityShape(mesh->getType(e));
+    parent = p;
+    nen = shape->countNodes();
+    nc = f->countComponents();
+    getNodeData();
+  }
+  void getNodeData()
+  {
+    reinterpret_cast<FieldDataOf<T>*>(field->getData())->getElementData(entity,nodeData);
+  }
+  NewArray<T> nodeData;
+  FieldBase* field;
+  Mesh* mesh;
+  MeshEntity* entity;
+  EntityShape* shape;
+  VectorElement* parent;
+  int nen;
+  int nc;
 };
 
-Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim);
+
+template <class T>
+int countNodes(ElementBase<T> * e)
+{
+  return e->getShape()->countNodes();
+}
+
+template <class T>
+void getShapeValues(ElementBase<T> * e, Vector3 const& local, NewArray<double>& values)
+{
+  e->getShape()->getValues(e->getMesh(),e->getEntity(),local,values);
+}
+
+template <class T>
+void getShapeGrads(ElementBase<T> * e, Vector3 const& local, NewArray<Vector3>& grads)
+{
+  e->getGlobalGradients(local,grads);
+}
 
 }//namespace apf
 
diff --git a/apf/apfElementOf.h b/apf/apfElementOf.h
index 9b20565d7..fe0aabc75 100644
--- a/apf/apfElementOf.h
+++ b/apf/apfElementOf.h
@@ -37,11 +37,11 @@ class ElementOf : public Element
       getComponents(local, reinterpret_cast<double*>(value));
       return value[0];
     }
-    void getValues(NewArray<T>& values)
+    void getValues(NewArray<T>& values, int nc = 1)
     {
-      values.allocate(nen);
+      values.allocate(nen * nc);
       T* nodeValues = getNodeValues();
-      for (int i=0; i < nen; ++i)
+      for (int i=0; i < nen * nc; ++i)
         values[i] = nodeValues[i];
     }
 };
diff --git a/apf/apfElementType.h b/apf/apfElementType.h
new file mode 100644
index 000000000..681a0d402
--- /dev/null
+++ b/apf/apfElementType.h
@@ -0,0 +1,16 @@
+#ifndef APFELEMENTTYPE_H
+#define APFELEMENTTYPE_H
+
+#include "apfComplexType.h"
+
+namespace apf
+{
+
+template <class T>
+class ElementBase;
+typedef ElementBase<double> Element;
+typedef ElementBase<double_complex> ComplexElement;
+
+};
+
+#endif
diff --git a/apf/apfField.cc b/apf/apfField.cc
index 2196482ca..babe2e8ae 100644
--- a/apf/apfField.cc
+++ b/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,35 +97,15 @@ 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);
 }
 
+bool isFrozen(FieldBase * fb)
+{
+  return fb->getData()->isFrozen();
+}
+
 } //namespace apf
diff --git a/apf/apfField.h b/apf/apfField.h
index c68fc50a8..6b5703871 100644
--- a/apf/apfField.h
+++ b/apf/apfField.h
@@ -10,10 +10,10 @@
 
 #include <string>
 #include "apfMesh.h"
+#include "apfElementType.h"
 
 namespace apf {
 
-class Element;
 class FieldData;
 
 class FieldBase
@@ -71,12 +71,14 @@ 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);
 
 } //namespace apf
 
diff --git a/apf/apfFieldData.cc b/apf/apfFieldData.cc
index 2cc1a2c8b..c85b80572 100644
--- a/apf/apfFieldData.cc
+++ b/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>;
@@ -294,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
diff --git a/apf/apfFieldOf.h b/apf/apfFieldOf.h
index 209c04555..765663f9e 100644
--- a/apf/apfFieldOf.h
+++ b/apf/apfFieldOf.h
@@ -13,6 +13,7 @@
 
 namespace apf {
 
+
 template <class T>
 class FieldOf;
 
@@ -23,6 +24,9 @@ void project(FieldOf<T>* to, FieldOf<T>* from);
 template <class T>
 void axpy(double a, FieldOf<T>* x, FieldOf<T>* y);
 
+// template <class T, typename = void>
+// class FieldOf;
+// typename std::enable_if<std::is_standard_layout<T>::value>::type
 template <class T>
 class FieldOf : public Field
 {
@@ -40,14 +44,13 @@ class FieldOf : public Field
     }
     void project(Field* from)
     {
-      apf::project<T>(this,static_cast<FieldOf<T>*>(from));
+      apf::project(this,static_cast<FieldOf<T>*>(from));
     }
     void axpy(double a, Field* x)
     {
-      apf::axpy<T>(a,static_cast<FieldOf<T>*>(x),this);
+      apf::axpy(a,static_cast<FieldOf<T>*>(x),this);
     }
 };
 
 }
-
 #endif
diff --git a/apf/apfMesh.cc b/apf/apfMesh.cc
index 02a34f193..2a7dd0dbd 100644
--- a/apf/apfMesh.cc
+++ b/apf/apfMesh.cc
@@ -8,6 +8,8 @@
 #include <PCU.h>
 #include "apfCoordData.h"
 #include "apfVectorField.h"
+#include "apfComplex.h"
+#include "apfComplexField.h"
 #include "apfShape.h"
 #include "apfNumbering.h"
 #include "apfTagData.h"
@@ -290,19 +292,40 @@ void Mesh::addField(Field* f)
   fields.push_back(f);
 }
 
+void Mesh::addComplexField(ComplexField* f)
+{
+  ifields.push_back(f);
+}
+
 void Mesh::removeField(Field* f)
 {
-  std::vector<Field*>::iterator it = std::find(fields.begin(),fields.end(),f);
+  auto it = std::find(fields.begin(),fields.end(),f);
   if (it != fields.end())
     fields.erase(it);
 }
 
+void Mesh::removeComplexField(ComplexField* f)
+{
+  auto it = std::find(ifields.begin(),ifields.end(),f);
+  if(it != ifields.end())
+    ifields.erase(it);
+}
+
 Field* Mesh::findField(const char* name)
 {
   std::string n(name);
   for (size_t i=0; i < fields.size(); ++i)
-    if (n==getName(fields[i]))
-      return fields[i];
+    if (n == fields[i]->getName())
+      return dynamic_cast<Field*>(fields[i]);
+  return 0;
+}
+
+ComplexField* Mesh::findComplexField(const char* name)
+{
+  std::string n(name);
+  for(auto fld : ifields)
+    if(n == fld->getName())
+      return fld;
   return 0;
 }
 
@@ -311,11 +334,21 @@ int Mesh::countFields()
   return static_cast<int>(fields.size());
 }
 
+int Mesh::countComplexFields()
+{
+  return static_cast<int>(ifields.size());
+}
+
 Field* Mesh::getField(int i)
 {
   return fields[i];
 }
 
+ComplexField* Mesh::getComplexField(int i)
+{
+  return ifields[i];
+}
+
 void Mesh::addNumbering(Numbering* n)
 {
   numberings.push_back(n);
@@ -727,13 +760,12 @@ void changeMeshShape(Mesh* m, FieldShape* newShape, bool project)
   m->changeShape(newShape, project);
 }
 
-void unfreezeFields(Mesh* m) {
-  Field* f;
-  for (int i=0; i<m->countFields(); i++) {
-    f = m->getField(i);
-    if (isFrozen(f)) 
-      unfreeze(f);
-  }
+void unfreezeFields(Mesh* m)
+{
+  for(int ii = 0; ii < m->countFields(); ++ii)
+    unfreeze(m->getField(ii));
+  for(int ii = 0; ii < m->countComplexFields(); ++ii)
+    unfreeze(m->getComplexField(ii));
   m->hasFrozenFields = false;
 }
 
diff --git a/apf/apfMesh.h b/apf/apfMesh.h
index bdba056a1..ae9714bc8 100644
--- a/apf/apfMesh.h
+++ b/apf/apfMesh.h
@@ -14,6 +14,7 @@
 #include <vector>
 #include <map>
 #include <set>
+#include "apfComplexType.h"
 #include "apfVector.h"
 #include "apfDynamicArray.h"
 
@@ -21,8 +22,10 @@ struct gmi_model;
 
 namespace apf {
 
+class ComplexField;
 class FieldShape;
 class Field;
+class FieldBase;
 template <class T>
 class NumberingOf;
 /** \brief Numbering is meant to be a 32-bit local numbering */
@@ -225,6 +228,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 +244,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 +272,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;
@@ -356,16 +367,21 @@ class Mesh
       \details most users don't need this, functions in apf.h
                automatically call it */
     void addField(Field* f);
+    void addComplexField(ComplexField* f);
     /** \brief disassociate a field from this mesh
       \details most users don't need this, functions in apf.h
                automatically call it */
     void removeField(Field* f);
+    void removeComplexField(ComplexField* f);
     /** \brief lookup a field by its unique name */
     Field* findField(const char* name);
+    ComplexField* findComplexField(const char* name);
     /** \brief get the number of associated fields */
     int countFields();
+    int countComplexFields();
     /** \brief get the i'th associated field */
     Field* getField(int i);
+    ComplexField* getComplexField(int i);
     /** \brief associate a numbering with this mesh
       \details most users don't need this, functions in apfNumbering.h
                automatically call it */
@@ -392,6 +408,7 @@ class Mesh
   protected:
     Field* coordinateField;
     std::vector<Field*> fields;
+    std::vector<ComplexField*> ifields;
     std::vector<Numbering*> numberings;
     std::vector<GlobalNumbering*> globalNumberings;
 };
diff --git a/apf/apfTagData.h b/apf/apfTagData.h
index 547be15a3..8b3b5bc1d 100644
--- a/apf/apfTagData.h
+++ b/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
 {
diff --git a/apf/apfZero.cc b/apf/apfZero.cc
new file mode 100644
index 000000000..c019cfbec
--- /dev/null
+++ b/apf/apfZero.cc
@@ -0,0 +1,18 @@
+#include "apfZero.h"
+#include "apfComplex.h"
+#include "apfComplexField.h"
+namespace apf
+{
+template <>
+void setComponents<double>(FieldBase* f, MeshEntity* e, int node, double const * components)
+{
+  setComponents(static_cast<Field*>(f),e,node,components);
+}
+
+template <>
+void setComponents<double_complex>(FieldBase* f, MeshEntity* e, int node, double_complex const * components)
+{
+  setComponents(static_cast<ComplexField*>(f),e,node,components);
+}
+
+}
diff --git a/apf/apfZero.h b/apf/apfZero.h
new file mode 100644
index 000000000..bfee2bfdc
--- /dev/null
+++ b/apf/apfZero.h
@@ -0,0 +1,40 @@
+#ifndef APFZERO_H
+#define APFZERO_H
+
+#include "apfField.h"
+#include "apfNew.h"
+
+namespace apf
+{
+
+template <class T>
+void setComponents(FieldBase* f, MeshEntity* e, int node, T const * components);
+
+template <class T>
+struct ZeroOp : public FieldOp
+{
+  ZeroOp(FieldBase * fb)
+  {
+    f = fb;
+    int cmps = fb->countComponents();
+    data.allocate(cmps);
+    for (int i = 0; i < cmps; ++i)
+      data[i] = 0;
+    ent = 0;
+  }
+  bool inEntity(MeshEntity* e)
+  {
+    ent = e;
+    return true;
+  }
+  void atNode(int node)
+  {
+    setComponents<T>(f,ent,node,&data[0]);
+  }
+  FieldBase* f;
+  MeshEntity* ent;
+  apf::NewArray<T> data;
+};
+}
+
+#endif
diff --git a/apf/pkg_tribits.cmake b/apf/pkg_tribits.cmake
index a0c6cff0c..e7de49943 100644
--- a/apf/pkg_tribits.cmake
+++ b/apf/pkg_tribits.cmake
@@ -49,6 +49,8 @@ set(APF_SOURCES
   apfBoundaryToElementXi.cc
   apfSimplexAngleCalcs.cc
   apfFile.cc
+  apfComplex.h
+  apfComplexType.h
 )
 
 set(APF_HEADERS
@@ -73,6 +75,7 @@ set(APF_HEADERS
   apfField.h
   apfFieldData.h
   apfNumberingClass.h
+  apfElementType.h
 )
 
 set(APF_SOURCES
diff --git a/crv/crvQuality.cc b/crv/crvQuality.cc
index 04119933a..4daefd74e 100644
--- a/crv/crvQuality.cc
+++ b/crv/crvQuality.cc
@@ -142,26 +142,26 @@ static double getTetPartialJacobianDet(apf::NewArray<apf::Vector3>& nodes,
  *
  */
 static double Nijk(apf::NewArray<apf::Vector3>& nodes,
-    int d, int I, int J)
+    int d, int II, int J)
 {
   double sum = 0.;
-  int CD = trinomial(2*(d-1),I,J);
+  int CD = trinomial(2*(d-1),II,J);
   for(int j1 = 0; j1 <= J; ++j1){
-    int i1start = std::max(0,I+J-j1-(d-1));
-    int i1end = std::min(I,d-1-j1);
+    int i1start = std::max(0,II+J-j1-(d-1));
+    int i1end = std::min(II,d-1-j1);
     for(int i1 = i1start; i1 <= i1end; ++i1){
-      sum += trinomial(d-1,i1,j1)*trinomial(d-1,I-i1,J-j1)
-          *getTriPartialJacobianDet(nodes,d,i1,j1,I-i1,J-j1);
+      sum += trinomial(d-1,i1,j1)*trinomial(d-1,II-i1,J-j1)
+          *getTriPartialJacobianDet(nodes,d,i1,j1,II-i1,J-j1);
     }
   }
   return sum*d*d/CD;
 }
 
 static double Nijkl(apf::NewArray<apf::Vector3>& nodes,
-    int d, int I, int J, int K)
+    int d, int II, int J, int K)
 {
   double sum = 0.;
-  int CD = quadnomial(3*(d-1),I,J,K);
+  int CD = quadnomial(3*(d-1),II,J,K);
 
   for(int k1 = 0; k1 <= K; ++k1){
     int k2start = std::max(0,K-k1-(d-1));
@@ -169,12 +169,12 @@ static double Nijkl(apf::NewArray<apf::Vector3>& nodes,
       for (int j1 = 0; j1 <= J; ++j1){
         int j2start = std::max(0,J-j1-(d-1));
         for (int j2 = j2start; j2 <= J-j1; ++j2){
-          int i1end = std::min(I,d-1-j1-k1);
+          int i1end = std::min(II,d-1-j1-k1);
           for (int i1 = 0; i1 <= i1end; ++i1){
-            int i2start = std::max(0,I+J+K-i1-j1-k1-j2-k2-(d-1));
-            int i2end = std::min(I-i1,d-1-j2-k2);
+            int i2start = std::max(0,II+J+K-i1-j1-k1-j2-k2-(d-1));
+            int i2end = std::min(II-i1,d-1-j2-k2);
             for (int i2 = i2start; i2 <= i2end; ++i2){
-              int i3 = I-i1-i2;
+              int i3 = II-i1-i2;
               int j3 = J-j1-j2;
               int k3 = K-k1-k2;
               sum += quadnomial(d-1,i1,j1,k1)*quadnomial(d-1,i2,j2,k2)
@@ -215,9 +215,9 @@ static double calcMaxJacDet(int n, apf::NewArray<double>& nodes)
 static void getTriJacDetNodes(int P, apf::NewArray<apf::Vector3>& elemNodes,
     apf::NewArray<double>& nodes)
 {
-  for (int I = 0; I <= 2*(P-1); ++I)
-    for (int J = 0; J <= 2*(P-1)-I; ++J)
-      nodes[getTriNodeIndex(2*(P-1),I,J)] = Nijk(elemNodes,P,I,J);
+  for (int II = 0; II <= 2*(P-1); ++II)
+    for (int J = 0; J <= 2*(P-1)-II; ++J)
+      nodes[getTriNodeIndex(2*(P-1),II,J)] = Nijk(elemNodes,P,II,J);
 }
 //static void getTetJacDetNodes(int P, apf::NewArray<apf::Vector3>& elemNodes,
 //    apf::NewArray<double>& nodes)
@@ -581,11 +581,11 @@ double computeTriJacobianDetFromBezierFormulation(apf::Mesh* m,
   apf::NewArray<apf::Vector3> nodes;
   apf::getVectorNodes(elem,nodes);
 
-  for (int I = 0; I <= 2*(P-1); ++I){
-    for (int J = 0; J <= 2*(P-1)-I; ++J){
-      detJ += trinomial(2*(P-1),I,J)
-          *Bijk(I,J,2*(P-1)-I-J,1.-xi[0]-xi[1],xi[0],xi[1])
-          *Nijk(nodes,P,I,J);
+  for (int II = 0; II <= 2*(P-1); ++II){
+    for (int J = 0; J <= 2*(P-1)-II; ++J){
+      detJ += trinomial(2*(P-1),II,J)
+          *Bijk(II,J,2*(P-1)-II-J,1.-xi[0]-xi[1],xi[0],xi[1])
+          *Nijk(nodes,P,II,J);
     }
   }
   apf::destroyElement(elem);
@@ -601,12 +601,12 @@ double computeTetJacobianDetFromBezierFormulation(apf::Mesh* m,
   apf::getVectorNodes(elem,nodes);
 
   double detJ = 0.;
-  for (int I = 0; I <= 3*(P-1); ++I){
-    for (int J = 0; J <= 3*(P-1)-I; ++J){
-      for (int K = 0; K <= 3*(P-1)-I-J; ++K){
-        detJ += quadnomial(3*(P-1),I,J,K)*Bijkl(I,J,K,3*(P-1)-I-J-K,
+  for (int II = 0; II <= 3*(P-1); ++II){
+    for (int J = 0; J <= 3*(P-1)-II; ++J){
+      for (int K = 0; K <= 3*(P-1)-II-J; ++K){
+        detJ += quadnomial(3*(P-1),II,J,K)*Bijkl(II,J,K,3*(P-1)-II-J-K,
             1.-xi[0]-xi[1]-xi[2],xi[0],xi[1],xi[2])
-            *Nijkl(nodes,P,I,J,K);
+            *Nijkl(nodes,P,II,J,K);
       }
     }
   }
diff --git a/mds/apfMDS.cc b/mds/apfMDS.cc
index 322322209..f1e9d3473 100644
--- a/mds/apfMDS.cc
+++ b/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);
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index e1562cef1..02829ba78 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -134,6 +134,7 @@ test_exe_func(integrate integrate.cc)
 test_exe_func(align align.cc)
 test_exe_func(field_io field_io.cc)
 test_exe_func(tensor tensor.cc)
+test_exe_func(complex complex.cc)
 test_exe_func(test_AD test_AD.cc)
 test_exe_func(spr_test spr_test.cc)
 test_exe_func(reposition reposition.cc)
diff --git a/test/complex.cc b/test/complex.cc
new file mode 100644
index 000000000..d4f9c6b07
--- /dev/null
+++ b/test/complex.cc
@@ -0,0 +1,100 @@
+#include <apf.h>
+#include <apfComplex.h>
+#include <apfMDS.h>
+#include <apfMesh2.h>
+#include <apfShape.h>
+#include <gmi_mesh.h>
+#include <PCU.h>
+#include <lionPrint.h>
+#include <pcu_util.h>
+
+int main(int ac, char** av)
+{
+  // requires C++ complex, not C complex
+  using namespace std::complex_literals;
+  PCU_ALWAYS_ASSERT(ac == 3);
+  MPI_Init(&ac,&av);
+  PCU_Comm_Init();
+  gmi_register_mesh();
+  apf::Mesh2* m = apf::loadMdsMesh(av[1], av[2]);
+  static const int cmps = 2;
+  apf::ComplexField* f = apf::createComplexPackedField(m,"foo",cmps);
+  int vrts = 0;
+  apf::zeroField(f);
+  double_complex zero = {0,0};
+  double_complex r[cmps];
+  double_complex w[cmps] = { {1,2},{2,1} };
+  {
+    apf::MeshEntity* vert = NULL;
+    apf::MeshIterator* it = m->begin(0);
+    while ((vert = m->iterate(it)))
+    {
+      apf::getComponents(f,vert,0,&r[0]);
+      PCU_ALWAYS_ASSERT(r[0] == zero && r[1] == zero);
+      apf::setComponents(f,vert,0,&w[0]);
+      apf::getComponents(f,vert,0,&r[0]);
+      PCU_ALWAYS_ASSERT(r[0] == w[0] && r[1] == w[1]);
+      ++vrts;
+    }
+    m->end(it);
+  }
+  apf::freeze(f);
+  PCU_ALWAYS_ASSERT(apf::isFrozen(f));
+  {
+    apf::MeshEntity* vert = NULL;
+    apf::MeshIterator* it = m->begin(0);
+    while ((vert = m->iterate(it)))
+    {
+      apf::getComponents(f,vert,0,&r[0]);
+      PCU_ALWAYS_ASSERT(r[0] == w[0] && r[1] == w[1]);
+    }
+    m->end(it);
+  }
+  double_complex * foo_arr = apf::getComplexArrayData(f);
+  for(int vv = 0; vv < vrts; ++vv)
+    PCU_ALWAYS_ASSERT(foo_arr[cmps*vv] == r[0] && foo_arr[cmps*vv+1] == r[1]);
+  apf::unfreeze(f);
+  foo_arr = NULL; // im certain if we tried to access this we would get an error,
+                  // but there isn't a good way to prevent this without changing the API
+  PCU_ALWAYS_ASSERT(!apf::isFrozen(f));
+
+  // TODO: test getComponents at a parametric location
+  //       getshapevalues and getshapegrads only really use the meshelement so dont bother
+  {
+    apf::MeshEntity* e = NULL;
+    int dim = m->getDimension();
+    apf::MeshIterator* it = m->begin(dim);
+    while ((e = m->iterate(it)))
+    {
+      apf::MeshElement * me = apf::createMeshElement(m,e);
+      apf::ComplexElement * el1 = apf::createElement(f,me);
+      apf::ComplexElement * el2 = apf::createElement(f,e);
+      apf::MeshElement * me1 = apf::getMeshElement(el1);
+      apf::MeshElement * me2 = apf::getMeshElement(el2);
+      PCU_ALWAYS_ASSERT(me1 == me && me2 == NULL);
+      apf::MeshEntity * e1 = apf::getMeshEntity(el1);
+      apf::MeshEntity * e2 = apf::getMeshEntity(el2);
+      PCU_ALWAYS_ASSERT(e1 == e && e2 == e);
+      int nds1 = apf::countNodes(el1);
+      int nds2 = apf::countNodes(el2);
+      PCU_ALWAYS_ASSERT(nds1 == nds2);
+      apf::NewArray<double_complex> dofs1(cmps*nds1);
+      apf::NewArray<double_complex> dofs2(cmps*nds2);
+      apf::getPackedNodes(el1,dofs1);
+      apf::getPackedNodes(el2,dofs2);
+      for(int nd = 0; nd < nds1; ++nd)
+      {
+        PCU_ALWAYS_ASSERT(dofs1[nd*cmps] == w[0] && dofs1[nd*cmps+1] == w[1]);
+        PCU_ALWAYS_ASSERT(dofs2[nd*cmps] == w[0] && dofs2[nd*cmps+1] == w[1]);
+      }
+      apf::destroyElement(el2);
+      apf::destroyElement(el1);
+      apf::destroyMeshElement(me);
+    }
+    m->end(it);
+  }
+  m->destroyNative();
+  apf::destroyMesh(m);
+  PCU_Comm_Free();
+  MPI_Finalize();
+}
diff --git a/test/quality.cc b/test/quality.cc
index d4b675cde..2b9edd1ee 100644
--- a/test/quality.cc
+++ b/test/quality.cc
@@ -54,10 +54,10 @@ static void processMesh(apf::Mesh2* m)
   int d = m->getDimension();
   apf::MeshEntity* elem;
   apf::MeshIterator* elems = m->begin(d);
-  ma::IdentitySizeField I(m);
+  ma::IdentitySizeField II(m);
   while ((elem = m->iterate(elems)))
   {
-    double q = ma::measureElementQuality(m, &I, elem);
+    double q = ma::measureElementQuality(m, &II, elem);
     q = pow(q, (1.0/d));
     avgQuality += q;
     minQuality = fmin(minQuality, q);
diff --git a/test/testing.cmake b/test/testing.cmake
index a619319f2..b3501d5b9 100644
--- a/test/testing.cmake
+++ b/test/testing.cmake
@@ -34,6 +34,10 @@ mpi_test(integrate 1 ./integrate)
 mpi_test(qr_test 1 ./qr)
 mpi_test(base64 1 ./base64)
 mpi_test(tensor_test 1 ./tensor)
+mpi_test(complex_test 1
+         ./complex
+         "${MESHES}/cube/cube.dmg"
+         "${MESHES}/cube/pumi11/cube.smb")
 mpi_test(verify_convert 1 ./verify_convert)
 mpi_test(test_integrator 1
          ./test_integrator