From e72dfe2f6eeada1d0bb99b5d9cb11975400aa633 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Fri, 28 Jun 2019 19:41:27 -0400 Subject: [PATCH 01/17] changing and adding underlying types, at a point where a decision needs to be made about apf::Field becuase it assumes that the data it is storing is of type double --- apf/apf.cc | 3 +++ apf/apf.h | 10 +++++++++- apf/apfArrayData.cc | 21 ++++++++++++++++----- apf/apfFieldData.cc | 3 +++ apf/apfMesh.h | 11 ++++++++++- apf/apfTagData.h | 18 ++++++++++++++++++ mds/apfMDS.cc | 15 +++++++++++++++ pumi/pumi_ghost.cc | 2 +- 8 files changed, 75 insertions(+), 8 deletions(-) diff --git a/apf/apf.cc b/apf/apf.cc index c7b60a1a2..2da16b5c7 100644 --- a/apf/apf.cc +++ b/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" @@ -73,6 +74,8 @@ Field* makeField( f = new MatrixField(); else if (valueType == PACKED) f = new PackedField(components); + else if (valueType == COMPLEX_PACKED) + f = new ComplexPackedField(components); else fail("invalid valueType in field construction\n"); f->init(name,m,shape,data); diff --git a/apf/apf.h b/apf/apf.h index 1ea32d2df..9f20e55e9 100644 --- a/apf/apf.h +++ b/apf/apf.h @@ -133,6 +133,8 @@ enum ValueType { MATRIX, /** \brief a user-defined set of components */ PACKED, + /** \brief a user-defined set of complex-components */ + COMPLEX_PACKED, /** \brief placeholder used to set array sizes */ VALUE_TYPES }; @@ -704,6 +706,9 @@ bool isPrintable(Field* f); */ void fail(const char* why) __attribute__((noreturn)); + + + /** \brief Convert a Field from Tag to array storage. */ void freeze(Field* f); @@ -713,12 +718,15 @@ void unfreeze(Field* f); /** \brief Returns true iff the Field uses array storage. */ bool isFrozen(Field* f); +template +T* getArrayDataT(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. */ -double* getArrayData(Field* f); +inline double * getArrayData(Field* f) { return getArrayDataT(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..8e9021e59 100644 --- a/apf/apfArrayData.cc +++ b/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(FieldBase* field); template void freezeFieldData(FieldBase* field); template void freezeFieldData(FieldBase* field); +template void unfreezeFieldData(FieldBase * field); template void unfreezeFieldData(FieldBase* field); template void unfreezeFieldData(FieldBase* field); -double* getArrayData(Field* f) { - if (!isFrozen(f)) { +template +T* getArrayDataT(FieldBase* f) +{ + if(!isFrozen(f)) return 0; - } else { - FieldDataOf* p = f->getData(); - ArrayDataOf* a = static_cast* > (p); + else + { + FieldDataOf* p = reinterpret_cast*>(f->getData()); + ArrayDataOf* a = static_cast*>(p); return a->getDataArray(); } } +template double_complex* getArrayDataT(FieldBase* field); +template int* getArrayDataT(FieldBase* field); +template double* getArrayDataT(FieldBase* field); + + } diff --git a/apf/apfFieldData.cc b/apf/apfFieldData.cc index 2cc1a2c8b..66ac656c3 100644 --- a/apf/apfFieldData.cc +++ b/apf/apfFieldData.cc @@ -1,4 +1,5 @@ #include +#include "apfComplex.h" #include "apfFieldData.h" #include "apfShape.h" #include @@ -72,6 +73,7 @@ void synchronizeFieldData(FieldDataOf* data, Sharing* shr, bool delete_shr) } /* instantiate here */ +template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); template void synchronizeFieldData(FieldDataOf*, Sharing*, bool); @@ -250,6 +252,7 @@ int FieldDataOf::getElementData(MeshEntity* entity, NewArray& data) return n; } +template class FieldDataOf; template class FieldDataOf; template class FieldDataOf; template class FieldDataOf; diff --git a/apf/apfMesh.h b/apf/apfMesh.h index bdba056a1..be6af9ead 100644 --- a/apf/apfMesh.h +++ b/apf/apfMesh.h @@ -14,6 +14,7 @@ #include #include #include +#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 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; 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 : public TagMaker } }; +template <> +class TagHelper : 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 : public TagMaker { 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(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(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/pumi/pumi_ghost.cc b/pumi/pumi_ghost.cc index af59bfd63..5ada109ed 100644 --- a/pumi/pumi_ghost.cc +++ b/pumi/pumi_ghost.cc @@ -400,7 +400,7 @@ void pumi_ghost_create(pMesh m, Ghosting* plan) std::vector frozen_fields; for (int i=0; icountFields(); ++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 From c8ac116dcb77ce884b8e3c362288990dccd5998a Mon Sep 17 00:00:00 2001 From: William Tobin Date: Mon, 1 Jul 2019 09:38:48 -0400 Subject: [PATCH 02/17] adding complex field files and accomodating complex field freezing/unfreezing --- apf/CMakeLists.txt | 2 ++ apf/apf.cc | 2 +- apf/apfComplex.h | 14 ++++++++++++++ apf/apfComplexField.cc | 21 +++++++++++++++++++++ apf/apfComplexField.h | 25 +++++++++++++++++++++++++ apf/apfField.cc | 5 +++++ apf/apfField.h | 2 ++ 7 files changed, 70 insertions(+), 1 deletion(-) create mode 100644 apf/apfComplex.h create mode 100644 apf/apfComplexField.cc create mode 100644 apf/apfComplexField.h diff --git a/apf/CMakeLists.txt b/apf/CMakeLists.txt index 881a9fe65..ecafd44e2 100644 --- a/apf/CMakeLists.txt +++ b/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 diff --git a/apf/apf.cc b/apf/apf.cc index 2da16b5c7..fef999f2c 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -476,7 +476,7 @@ void unfreeze(Field* f) bool isFrozen(Field* f) { - return f->getData()->isFrozen(); + return isFrozen(static_cast(f)); } Function::~Function() diff --git a/apf/apfComplex.h b/apf/apfComplex.h new file mode 100644 index 000000000..640561b16 --- /dev/null +++ b/apf/apfComplex.h @@ -0,0 +1,14 @@ +#ifndef APFCOMPLEX_H_ +#define APFCOMPLEX_H_ + +#ifdef C_COMPLEX + #include +#endif + +#define CXX_COMPLEX 1 +#ifdef CXX_COMPLEX + #include + using double_complex = std::complex; +#endif + +#endif diff --git a/apf/apfComplexField.cc b/apf/apfComplexField.cc new file mode 100644 index 000000000..c62419601 --- /dev/null +++ b/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"); +} + +} diff --git a/apf/apfComplexField.h b/apf/apfComplexField.h new file mode 100644 index 000000000..298118043 --- /dev/null +++ b/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 diff --git a/apf/apfField.cc b/apf/apfField.cc index 2196482ca..f6b1d5a10 100644 --- a/apf/apfField.cc +++ b/apf/apfField.cc @@ -126,4 +126,9 @@ void zeroField(Field* 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..62bea2a1b 100644 --- a/apf/apfField.h +++ b/apf/apfField.h @@ -78,6 +78,8 @@ Field* makeField( FieldShape* shape, FieldData* data); +bool isFrozen(FieldBase * fb); + } //namespace apf #endif From dc9324c0ebafe759c9ed138c0ebdabe447225d42 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Mon, 8 Jul 2019 08:08:16 -0400 Subject: [PATCH 03/17] remove template function from API, add some type-checking to various field get/set functions --- apf/apf.cc | 12 ++++++++++-- apf/apf.h | 29 +++++++++++++++++++++-------- 2 files changed, 31 insertions(+), 10 deletions(-) diff --git a/apf/apf.cc b/apf/apf.cc index fef999f2c..ecd57d25f 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -74,8 +74,6 @@ Field* makeField( f = new MatrixField(); else if (valueType == PACKED) f = new PackedField(components); - else if (valueType == COMPLEX_PACKED) - f = new ComplexPackedField(components); else fail("invalid valueType in field construction\n"); f->init(name,m,shape,data); @@ -163,6 +161,7 @@ void destroyField(Field* f) void setScalar(Field* f, MeshEntity* e, int node, double value) { + PCU_DEBUG_ASSERT(f->getValueType() == SCALAR); ScalarField* field = static_cast(f); double tmp[1] = {value}; field->setNodeValue(e,node,tmp); @@ -170,6 +169,7 @@ void setScalar(Field* f, MeshEntity* e, int node, double value) double getScalar(Field* f, MeshEntity* e, int node) { + PCU_DEBUG_ASSERT(f->getValueType() == SCALAR); ScalarField* field = static_cast(f); double value[1]; field->getNodeValue(e,node,value); @@ -178,6 +178,7 @@ 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(f); Vector3 tmp[1] = {value}; field->setNodeValue(e,node,tmp); @@ -185,6 +186,7 @@ void setVector(Field* f, MeshEntity* e, int node, Vector3 const& value) void getVector(Field* f, MeshEntity* e, int node, Vector3& value) { + PCU_DEBUG_ASSERT(f->getValueType() == VECTOR); VectorField* field = static_cast(f); Vector3 tmp[1]; field->getNodeValue(e,node,tmp); @@ -193,6 +195,7 @@ 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(f); Matrix3x3 tmp[1] = {value}; field->setNodeValue(e,node,tmp); @@ -200,6 +203,7 @@ void setMatrix(Field* f, MeshEntity* e, int node, Matrix3x3 const& value) void getMatrix(Field* f, MeshEntity* e, int node, Matrix3x3& value) { + PCU_DEBUG_ASSERT(f->getValueType() == MATRIX); MatrixField* field = static_cast(f); Matrix3x3 tmp[1]; field->getNodeValue(e,node,tmp); @@ -208,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); } diff --git a/apf/apf.h b/apf/apf.h index 9f20e55e9..28b94242b 100644 --- a/apf/apf.h +++ b/apf/apf.h @@ -11,6 +11,7 @@ #include "apfMatrix.h" #include "apfNew.h" #include "apfDynamicArray.h" +#include "apfComplex.h" #include #include @@ -133,8 +134,6 @@ enum ValueType { MATRIX, /** \brief a user-defined set of components */ PACKED, - /** \brief a user-defined set of complex-components */ - COMPLEX_PACKED, /** \brief placeholder used to set array sizes */ VALUE_TYPES }; @@ -706,9 +705,6 @@ bool isPrintable(Field* f); */ void fail(const char* why) __attribute__((noreturn)); - - - /** \brief Convert a Field from Tag to array storage. */ void freeze(Field* f); @@ -718,15 +714,32 @@ void unfreeze(Field* f); /** \brief Returns true iff the Field uses array storage. */ bool isFrozen(Field* f); -template -T* getArrayDataT(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, + 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. */ -inline double * getArrayData(Field* f) { return getArrayDataT(f); } +double_complex* getComplexArrayData(Field * f); /** \brief Initialize all nodal values with all-zero components */ void zeroField(Field* f); From 1f0910538c3775fec5557812614687158b7a00c9 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Mon, 8 Jul 2019 08:10:03 -0400 Subject: [PATCH 04/17] closing parenthesis, gotta squash this one for sure --- apf/apf.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/apf/apf.cc b/apf/apf.cc index ecd57d25f..16af91da7 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -213,14 +213,14 @@ 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->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->getValueType() == SCALAR || f->getValueType() == PACKED)); f->getData()->getNodeComponents(e,node,components); } From 246a8438efc5c526e95bd895d0eea8e56ca52a53 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Mon, 8 Jul 2019 12:38:54 -0400 Subject: [PATCH 05/17] enforcing a small amount of type safety to get/set calls in the API, implementing ability to get Element nodes for a packed field, restricting ElementOf to be standard layout so the underlying data retrieval cast is valid (see issue #239) --- apf/apf.cc | 32 ++++++++++++++++++++++++++++---- apf/apf.h | 16 +++++----------- apf/apfArrayData.cc | 16 ++++++++++++++++ apf/apfElementOf.h | 9 +++++---- 4 files changed, 54 insertions(+), 19 deletions(-) diff --git a/apf/apf.cc b/apf/apf.cc index 16af91da7..69e4cd30f 100644 --- a/apf/apf.cc +++ b/apf/apf.cc @@ -75,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; @@ -127,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) @@ -251,42 +255,49 @@ MeshEntity* getMeshEntity(Element* e) double getScalar(Element* e, Vector3 const& param) { + PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); ScalarElement* element = static_cast(e); return element->getValue(param); } void getGrad(Element* e, Vector3 const& param, Vector3& g) { + PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); ScalarElement* element = static_cast(e); element->grad(param,g); } void getVector(Element* e, Vector3 const& param, Vector3& value) { + PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); VectorElement* element = static_cast(e); value = element->getValue(param); } double getDiv(Element* e, Vector3 const& param) { + PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); VectorElement* element = static_cast(e); return element->div(param); } void getCurl(Element* e, Vector3 const& param, Vector3& c) { + PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); VectorElement* element = static_cast(e); return element->curl(param,c); } void getVectorGrad(Element* e, Vector3 const& param, Matrix3x3& deriv) { + PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); VectorElement* element = static_cast(e); return element->grad(param,deriv); } void getMatrix(Element* e, Vector3 const& param, Matrix3x3& value) { + PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); MatrixElement* element = static_cast(e); value = element->getValue(param); } @@ -294,6 +305,7 @@ void getMatrix(Element* e, Vector3 const& param, Matrix3x3& value) void getMatrixGrad(Element* e, Vector3 const& param, Vector<27>& deriv) { + PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); MatrixElement* element = static_cast(e); return element->grad(param,deriv); } @@ -382,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(e)); } void getScalarNodes(Element* e, NewArray& values) { + PCU_DEBUG_ASSERT(dynamic_cast*>(e) != NULL); ElementOf* element = static_cast*>(e); element->getValues(values); } void getVectorNodes(Element* e, NewArray& values) { + PCU_DEBUG_ASSERT(dynamic_cast*>(e) != NULL); ElementOf* element = static_cast*>(e); element->getValues(values); } void getMatrixNodes(Element* e, NewArray& values) { + PCU_DEBUG_ASSERT(dynamic_cast*>(e) != NULL); ElementOf* element = static_cast*>(e); element->getValues(values); } +void getPackedNodes(Element* e, NewArray& values) +{ + FieldBase* f = e->getFieldBase(); + int cmps = f->countComponents(); + PCU_DEBUG_ASSERT(dynamic_cast*>(e) != NULL); + ElementOf* element = static_cast*>(e); + element->getValues(values,cmps); +} + void getShapeValues(Element* e, Vector3 const& local, NewArray& values) { diff --git a/apf/apf.h b/apf/apf.h index 28b94242b..b517ac637 100644 --- a/apf/apf.h +++ b/apf/apf.h @@ -11,7 +11,6 @@ #include "apfMatrix.h" #include "apfNew.h" #include "apfDynamicArray.h" -#include "apfComplex.h" #include #include @@ -44,7 +43,7 @@ template class ReductionOp; template 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 @@ -570,6 +569,10 @@ void getVectorNodes(Element* e, NewArray& values); */ void getMatrixNodes(Element* e, NewArray& values); +/** \brief Returns the element nodal values for a packed field + */ +void getPackedNodes(Element* e, NewArray& values); + /** \brief Returns the shape function values at a point */ void getShapeValues(Element* e, Vector3 const& local, @@ -732,15 +735,6 @@ double* getArrayData(Field* f); */ 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); diff --git a/apf/apfArrayData.cc b/apf/apfArrayData.cc index 8e9021e59..dff22c60a 100644 --- a/apf/apfArrayData.cc +++ b/apf/apfArrayData.cc @@ -2,6 +2,8 @@ #include "apfComplex.h" #include "apfNumbering.h" #include "apfTagData.h" +#include +#include namespace apf { @@ -136,6 +138,17 @@ template void unfreezeFieldData(FieldBase* field); template 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::value) || + (scalar == Mesh::INT && std::is_same::value) || + (scalar == Mesh::LONG && std::is_same::value) || + (scalar == Mesh::COMPLEX && std::is_same::value) + ); if(!isFrozen(f)) return 0; else @@ -150,5 +163,8 @@ template double_complex* getArrayDataT(FieldBase* field); template int* getArrayDataT(FieldBase* field); template double* getArrayDataT(FieldBase* field); +double * getArrayData(Field * f) { return getArrayDataT(f); } +int * getIntArrayData(Field * f) { return getArrayDataT(f); } +double_complex * getComplexArrayData(Field * f) { return getArrayDataT(f); } } diff --git a/apf/apfElementOf.h b/apf/apfElementOf.h index 9b20565d7..104826f94 100644 --- a/apf/apfElementOf.h +++ b/apf/apfElementOf.h @@ -14,7 +14,8 @@ namespace apf { -template +template ::value, bool> = 0 > class ElementOf : public Element { public: @@ -37,11 +38,11 @@ class ElementOf : public Element getComponents(local, reinterpret_cast(value)); return value[0]; } - void getValues(NewArray& values) + void getValues(NewArray& 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]; } }; From 80c1c680a8904a5a54c9fc564e17b677c37a0156 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Mon, 8 Jul 2019 12:40:49 -0400 Subject: [PATCH 06/17] adding an ElementBase class and an intermediate ElementT class similar to the FieldBase->(Real)Field/ComplexField hierarchy --- apf/apfElement.cc | 35 ++++++------------ apf/apfElement.h | 90 ++++++++++++++++++++++++++++++++++------------- 2 files changed, 75 insertions(+), 50 deletions(-) diff --git a/apf/apfElement.cc b/apf/apfElement.cc index 82c89c602..193ca5f23 100644 --- a/apf/apfElement.cc +++ b/apf/apfElement.cc @@ -12,7 +12,7 @@ namespace apf { -void Element::init(Field* f, MeshEntity* e, VectorElement* p) +void ElementBase::init(FieldBase* f, MeshEntity* e, VectorElement* p) { field = f; mesh = f->getMesh(); @@ -24,20 +24,16 @@ void Element::init(Field* f, MeshEntity* e, VectorElement* p) getNodeData(); } -Element::Element(Field* f, MeshEntity* e) +ElementBase::ElementBase(FieldBase* f, MeshEntity* e) { init(f,e,0); } -Element::Element(Field* f, VectorElement* p) +ElementBase::ElementBase(FieldBase* f, VectorElement* p) { init(f,p->getEntity(),p); } -Element::~Element() -{ -} - Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim) { switch (dim) { @@ -72,8 +68,13 @@ Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim) } } -void Element::getGlobalGradients(Vector3 const& local, - NewArray& globalGradients) +int countNodes(ElementBase * e) +{ + return e->getShape()->countNodes(); +} + +void ElementBase::getGlobalGradients(Vector3 const& local, + NewArray& globalGradients) { Matrix3x3 J; parent->getJacobian(local,J); @@ -85,20 +86,4 @@ void Element::getGlobalGradients(Vector3 const& local, globalGradients[i] = jinv * localGradients[i]; } -void Element::getComponents(Vector3 const& xi, double* c) -{ - NewArray 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..797f91003 100644 --- a/apf/apfElement.h +++ b/apf/apfElement.h @@ -10,6 +10,7 @@ #include "apfMesh.h" #include "apfField.h" +#include "apfFieldData.h" #include "apfShape.h" namespace apf { @@ -17,36 +18,75 @@ namespace apf { class EntityShape; class VectorElement; -class Element +class ElementBase +{ +public: + ElementBase(FieldBase* f, MeshEntity* e); + ElementBase(FieldBase* f, VectorElement* p); + virtual ~ElementBase() {} + void getGlobalGradients(Vector3 const& local, + NewArray& 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; } +protected: + void init(FieldBase* f, MeshEntity* e, VectorElement* p); + virtual void getNodeData() = 0; + FieldBase* field; + Mesh* mesh; + MeshEntity* entity; + EntityShape* shape; + VectorElement* parent; + int nen; + int nc; +}; + +// not to be confused with ElementOf, which is about +// Scalar/Vector/Matrix node/value types +// this is about the underlying scalar type (real/complex) +template +class ElementT : public ElementBase +{ +public: + ElementT(FieldBase* f, MeshEntity* e) : ElementBase(f,e) { } + ElementT(FieldBase* f, VectorElement* p) : ElementBase(f,p) { } + virtual ~ElementT() {} + void getComponents(Vector3 const& xi, T * c) + { + NewArray 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: + virtual void getNodeData() + { + reinterpret_cast*>(field->getData())->getElementData(entity,nodeData); + } + NewArray nodeData; +}; + +class Element : public ElementT { public: - Element(Field* f, MeshEntity* e); - Element(Field* f, VectorElement* p); - virtual ~Element(); - void getGlobalGradients(Vector3 const& local, - NewArray& 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 nodeData; + Element(FieldBase* f, MeshEntity* e) + : ElementT(f,e) + { } + Element(FieldBase* f, VectorElement* p) + : ElementT(f,p) + { } }; Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim); +int countNodes(ElementBase * e); }//namespace apf From 07965d0ab517951f91829edf39bdcbf8b7ae4691 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Mon, 8 Jul 2019 12:42:34 -0400 Subject: [PATCH 07/17] implementation of the initial complex fields API, including extension of the Mesh API (though no downstream changes should be required) --- apf/CMakeLists.txt | 1 + apf/apfComplex.h | 53 +++++++++++++++++++ apf/apfComplexField.cc | 113 ++++++++++++++++++++++++++++++++++++++--- apf/apfComplexField.h | 36 ++++++++++--- apf/apfField.cc | 33 ++---------- apf/apfField.h | 12 ++--- apf/apfFieldData.cc | 2 + apf/apfFieldOf.h | 8 +-- apf/apfMesh.cc | 51 +++++++++++++++---- apf/apfMesh.h | 7 +++ apf/apfZero.cc | 18 +++++++ apf/apfZero.h | 39 ++++++++++++++ 12 files changed, 314 insertions(+), 59 deletions(-) create mode 100644 apf/apfZero.cc create mode 100644 apf/apfZero.h diff --git a/apf/CMakeLists.txt b/apf/CMakeLists.txt index ecafd44e2..0c853935e 100644 --- a/apf/CMakeLists.txt +++ b/apf/CMakeLists.txt @@ -48,6 +48,7 @@ set(SOURCES apfSimplexAngleCalcs.cc apfFile.cc apfMIS.cc + apfZero.cc ) # Package headers diff --git a/apf/apfComplex.h b/apf/apfComplex.h index 640561b16..ac00f2468 100644 --- a/apf/apfComplex.h +++ b/apf/apfComplex.h @@ -11,4 +11,57 @@ using double_complex = std::complex; #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 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& values); + +} + #endif diff --git a/apf/apfComplexField.cc b/apf/apfComplexField.cc index c62419601..47d810c2e 100644 --- a/apf/apfComplexField.cc +++ b/apf/apfComplexField.cc @@ -1,21 +1,122 @@ #include "apfComplexField.h" +#include "apfArrayData.h" #include "apfElement.h" +#include "apfTagData.h" +#include "apfZero.h" +#include -namespace apf { +namespace apf +{ + +FieldDataOf* ComplexField::getData() +{ + return static_cast*>(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); } +void freeze(ComplexField* f) +{ + if(isFrozen(f)) return; + f->getMesh()->hasFrozenFields = true; + freezeFieldData(f); +} + +void unfreeze(ComplexField* f) +{ + if(isFrozen(f)) + unfreezeFieldData(f); +} + +bool isFrozen(ComplexField* f) +{ + return isFrozen(static_cast(f)); +} + +void zero(ComplexField* f) +{ + ZeroOp 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(e)); +} + + + } diff --git a/apf/apfComplexField.h b/apf/apfComplexField.h index 298118043..bee7ed81f 100644 --- a/apf/apfComplexField.h +++ b/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 +{ +public: + ComplexElement(FieldBase* f, MeshEntity* e) + : ElementT(f,e) + { } + ComplexElement(FieldBase* f, VectorElement* p) + : ElementT(f,p) + { } +}; + +class ComplexField : public FieldBase +{ +public: + virtual ComplexElement * getElement(VectorElement*) = 0; + virtual int getValueType() const = 0; + virtual int getScalarType() { return Mesh::COMPLEX; } + FieldDataOf* 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; }; diff --git a/apf/apfField.cc b/apf/apfField.cc index f6b1d5a10..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,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 data; -}; - void zeroField(Field* f) { - ZeroOp op(f); + ZeroOp op(f); op.apply(f); } diff --git a/apf/apfField.h b/apf/apfField.h index 62bea2a1b..3659b44ce 100644 --- a/apf/apfField.h +++ b/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); diff --git a/apf/apfFieldData.cc b/apf/apfFieldData.cc index 66ac656c3..c85b80572 100644 --- a/apf/apfFieldData.cc +++ b/apf/apfFieldData.cc @@ -297,6 +297,8 @@ template void copyFieldData( FieldDataOf* from, FieldDataOf* to); template void copyFieldData( FieldDataOf* from, FieldDataOf* to); +template void copyFieldData( + FieldDataOf * from, FieldDataOf* to); template class MultiplyOp : public FieldOp diff --git a/apf/apfFieldOf.h b/apf/apfFieldOf.h index 209c04555..94bad5d7b 100644 --- a/apf/apfFieldOf.h +++ b/apf/apfFieldOf.h @@ -13,7 +13,10 @@ namespace apf { -template +template +using enable = typename std::enable_if_t::value, bool>; + +template class FieldOf; template @@ -24,7 +27,7 @@ template void axpy(double a, FieldOf* x, FieldOf* y); template -class FieldOf : public Field +class FieldOf::value>::type > : public Field { public: virtual ~FieldOf() {} @@ -49,5 +52,4 @@ class FieldOf : public Field }; } - #endif diff --git a/apf/apfMesh.cc b/apf/apfMesh.cc index 02a34f193..912d92829 100644 --- a/apf/apfMesh.cc +++ b/apf/apfMesh.cc @@ -8,6 +8,7 @@ #include #include "apfCoordData.h" #include "apfVectorField.h" +#include "apfComplexField.h" #include "apfShape.h" #include "apfNumbering.h" #include "apfTagData.h" @@ -290,19 +291,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::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(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 +333,21 @@ int Mesh::countFields() return static_cast(fields.size()); } +int Mesh::countComplexFields() +{ + return static_cast(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 +759,12 @@ void changeMeshShape(Mesh* m, FieldShape* newShape, bool project) m->changeShape(newShape, project); } -void unfreezeFields(Mesh* m) { - Field* f; - for (int i=0; icountFields(); 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 be6af9ead..178b32932 100644 --- a/apf/apfMesh.h +++ b/apf/apfMesh.h @@ -24,6 +24,7 @@ namespace apf { class FieldShape; class Field; +class FieldBase; template class NumberingOf; /** \brief Numbering is meant to be a 32-bit local numbering */ @@ -365,16 +366,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 */ @@ -401,6 +407,7 @@ class Mesh protected: Field* coordinateField; std::vector fields; + std::vector ifields; std::vector numberings; std::vector globalNumberings; }; 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(FieldBase* f, MeshEntity* e, int node, double const * components) +{ + setComponents(static_cast(f),e,node,components); +} + +template <> +void setComponents(FieldBase* f, MeshEntity* e, int node, double_complex const * components) +{ + setComponents(static_cast(f),e,node,components); +} + +} diff --git a/apf/apfZero.h b/apf/apfZero.h new file mode 100644 index 000000000..250ede35b --- /dev/null +++ b/apf/apfZero.h @@ -0,0 +1,39 @@ +#ifndef APFZERO_H +#define APFZERO_H + +#include "apfField.h" + +namespace apf +{ + +template +void setComponents(FieldBase* f, MeshEntity* e, int node, T const * components); + +template +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(f,ent,node,&data[0]); + } + FieldBase* f; + MeshEntity* ent; + apf::NewArray data; +}; +} + +#endif From 97a64129f765d63f1f3cd0d9f5687671c994f438 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Tue, 9 Jul 2019 18:17:16 -0400 Subject: [PATCH 08/17] quite a bit of work, reverting the type-checks on the API get/set funcitons, changing field zeroing and the element/elementof class hierarchy, adding complex field api to the mesh class, etc --- apf/apf.cc | 29 ++----------- apf/apf.h | 6 ++- apf/apfArrayData.cc | 7 ++- apf/apfComplex.h | 26 ++++------- apf/apfComplexField.cc | 53 +++++++++++++++------- apf/apfComplexField.h | 75 +++++++++++++++++++++++++++----- apf/apfComplexType.h | 10 +++++ apf/apfElement.cc | 51 +++++++++------------- apf/apfElement.h | 72 +++++++++++++++--------------- apf/apfElementOf.h | 5 ++- apf/apfField.h | 2 +- apf/apfFieldOf.h | 13 +++--- apf/apfMesh.cc | 1 + apf/apfMesh.h | 3 +- apf/apfZero.h | 1 + test/CMakeLists.txt | 1 + test/complex.cc | 99 ++++++++++++++++++++++++++++++++++++++++++ test/testing.cmake | 4 ++ 18 files changed, 310 insertions(+), 148 deletions(-) create mode 100644 apf/apfComplexType.h create mode 100644 test/complex.cc diff --git a/apf/apf.cc b/apf/apf.cc index 69e4cd30f..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" @@ -165,7 +166,6 @@ void destroyField(Field* f) void setScalar(Field* f, MeshEntity* e, int node, double value) { - PCU_DEBUG_ASSERT(f->getValueType() == SCALAR); ScalarField* field = static_cast(f); double tmp[1] = {value}; field->setNodeValue(e,node,tmp); @@ -173,7 +173,6 @@ void setScalar(Field* f, MeshEntity* e, int node, double value) double getScalar(Field* f, MeshEntity* e, int node) { - PCU_DEBUG_ASSERT(f->getValueType() == SCALAR); ScalarField* field = static_cast(f); double value[1]; field->getNodeValue(e,node,value); @@ -182,7 +181,6 @@ 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(f); Vector3 tmp[1] = {value}; field->setNodeValue(e,node,tmp); @@ -190,7 +188,6 @@ void setVector(Field* f, MeshEntity* e, int node, Vector3 const& value) void getVector(Field* f, MeshEntity* e, int node, Vector3& value) { - PCU_DEBUG_ASSERT(f->getValueType() == VECTOR); VectorField* field = static_cast(f); Vector3 tmp[1]; field->getNodeValue(e,node,tmp); @@ -199,7 +196,6 @@ 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(f); Matrix3x3 tmp[1] = {value}; field->setNodeValue(e,node,tmp); @@ -207,7 +203,6 @@ void setMatrix(Field* f, MeshEntity* e, int node, Matrix3x3 const& value) void getMatrix(Field* f, MeshEntity* e, int node, Matrix3x3& value) { - PCU_DEBUG_ASSERT(f->getValueType() == MATRIX); MatrixField* field = static_cast(f); Matrix3x3 tmp[1]; field->getNodeValue(e,node,tmp); @@ -216,15 +211,11 @@ 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); } @@ -255,49 +246,42 @@ MeshEntity* getMeshEntity(Element* e) double getScalar(Element* e, Vector3 const& param) { - PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); ScalarElement* element = static_cast(e); return element->getValue(param); } void getGrad(Element* e, Vector3 const& param, Vector3& g) { - PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); ScalarElement* element = static_cast(e); element->grad(param,g); } void getVector(Element* e, Vector3 const& param, Vector3& value) { - PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); VectorElement* element = static_cast(e); value = element->getValue(param); } double getDiv(Element* e, Vector3 const& param) { - PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); VectorElement* element = static_cast(e); return element->div(param); } void getCurl(Element* e, Vector3 const& param, Vector3& c) { - PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); VectorElement* element = static_cast(e); return element->curl(param,c); } void getVectorGrad(Element* e, Vector3 const& param, Matrix3x3& deriv) { - PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); VectorElement* element = static_cast(e); return element->grad(param,deriv); } void getMatrix(Element* e, Vector3 const& param, Matrix3x3& value) { - PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); MatrixElement* element = static_cast(e); value = element->getValue(param); } @@ -305,7 +289,6 @@ void getMatrix(Element* e, Vector3 const& param, Matrix3x3& value) void getMatrixGrad(Element* e, Vector3 const& param, Vector<27>& deriv) { - PCU_DEBUG_ASSERT(dynamic_cast(e) != NULL); MatrixElement* element = static_cast(e); return element->grad(param,deriv); } @@ -394,26 +377,23 @@ double computeCosAngle(Mesh* m, MeshEntity* pe, MeshEntity* e1, MeshEntity* e2, int countNodes(Element* e) { - return countNodes(static_cast(e)); + return countNodes(e); } void getScalarNodes(Element* e, NewArray& values) { - PCU_DEBUG_ASSERT(dynamic_cast*>(e) != NULL); ElementOf* element = static_cast*>(e); element->getValues(values); } void getVectorNodes(Element* e, NewArray& values) { - PCU_DEBUG_ASSERT(dynamic_cast*>(e) != NULL); ElementOf* element = static_cast*>(e); element->getValues(values); } void getMatrixNodes(Element* e, NewArray& values) { - PCU_DEBUG_ASSERT(dynamic_cast*>(e) != NULL); ElementOf* element = static_cast*>(e); element->getValues(values); } @@ -422,7 +402,6 @@ void getPackedNodes(Element* e, NewArray& values) { FieldBase* f = e->getFieldBase(); int cmps = f->countComponents(); - PCU_DEBUG_ASSERT(dynamic_cast*>(e) != NULL); ElementOf* element = static_cast*>(e); element->getValues(values,cmps); } @@ -430,13 +409,13 @@ void getPackedNodes(Element* e, NewArray& values) void getShapeValues(Element* e, Vector3 const& local, NewArray& values) { - e->getShape()->getValues(e->getMesh(), e->getEntity(), local,values); + getShapeValues(e,local,values); } void getShapeGrads(Element* e, Vector3 const& local, NewArray& grads) { - e->getGlobalGradients(local,grads); + getShapeGrads(e,local,grads); } FieldShape* getShape(Field* f) diff --git a/apf/apf.h b/apf/apf.h index b517ac637..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; @@ -724,7 +724,9 @@ bool isFrozen(Field* f); \note If the underlying field data type is NOT double, this will cause an assert fail in all compile modes. */ -double* getArrayData(Field* f); +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 diff --git a/apf/apfArrayData.cc b/apf/apfArrayData.cc index dff22c60a..0cf365542 100644 --- a/apf/apfArrayData.cc +++ b/apf/apfArrayData.cc @@ -163,8 +163,11 @@ template double_complex* getArrayDataT(FieldBase* field); template int* getArrayDataT(FieldBase* field); template double* getArrayDataT(FieldBase* field); -double * getArrayData(Field * f) { return getArrayDataT(f); } +double * getDoubleArrayData(Field * f) { return getArrayDataT(f); } int * getIntArrayData(Field * f) { return getArrayDataT(f); } -double_complex * getComplexArrayData(Field * f) { return getArrayDataT(f); } + +class ComplexField; +double_complex * getComplexArrayData(ComplexField * f) { return getArrayDataT(reinterpret_cast(f)); } + } diff --git a/apf/apfComplex.h b/apf/apfComplex.h index ac00f2468..39cf9a06a 100644 --- a/apf/apfComplex.h +++ b/apf/apfComplex.h @@ -1,22 +1,13 @@ #ifndef APFCOMPLEX_H_ #define APFCOMPLEX_H_ -#ifdef C_COMPLEX - #include -#endif - -#define CXX_COMPLEX 1 -#ifdef CXX_COMPLEX - #include - using double_complex = std::complex; -#endif +#include "apfComplexType.h" +#include "apfComplexElement.h" namespace apf { // forward decls for the interface -class ComplexElement; -class ComplexField; class Mesh; class FieldShape; class MeshEntity; @@ -26,11 +17,10 @@ class Vector3; template class NewArray; -ComplexField* createComplexField(Mesh* m, - const char* name, - int valueType, - int components, - FieldShape* shape); +ComplexField* createComplexPackedField(Mesh* m, + const char* name, + int components, + FieldShape* shape = NULL); void freeze(ComplexField* f); void unfreeze(ComplexField* f); @@ -58,8 +48,8 @@ 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 getShapeValues(ComplexElement* e, Vector3 const& local, NewArray& values); +void getShapeGrads(ComplexElement* e, Vector3 const& local, NewArray& grades); void getPackedNodes(ComplexElement* e, NewArray& values); } diff --git a/apf/apfComplexField.cc b/apf/apfComplexField.cc index 47d810c2e..457ee5d7d 100644 --- a/apf/apfComplexField.cc +++ b/apf/apfComplexField.cc @@ -1,6 +1,6 @@ #include "apfComplexField.h" #include "apfArrayData.h" -#include "apfElement.h" +#include "apfElementOf.h" #include "apfTagData.h" #include "apfZero.h" #include @@ -28,12 +28,12 @@ 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) +ComplexField* makeComplexGeneralField(Mesh* m, + const char* name, + int valueType, + int components, + FieldShape * shape, + FieldData * data) { PCU_ALWAYS_ASSERT(!m->findComplexField(name)); ComplexField* f = 0; @@ -48,12 +48,13 @@ ComplexField* makeComplexField(Mesh* m, ComplexField* createComplexPackedField(Mesh* m, const char* name, - int valueType, int components, FieldShape* shape) { - return makeComplexField(m,name,valueType,components,shape, - new TagDataOf); + if(shape == NULL) + shape = m->getShape(); + return makeComplexGeneralField(m,name,PACKED,components,shape, + new TagDataOf); } void freeze(ComplexField* f) @@ -74,7 +75,7 @@ bool isFrozen(ComplexField* f) return isFrozen(static_cast(f)); } -void zero(ComplexField* f) +void zeroField(ComplexField* f) { ZeroOp op(f); op.apply(f); @@ -82,13 +83,11 @@ void zero(ComplexField* 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); } @@ -107,6 +106,16 @@ 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); @@ -114,9 +123,23 @@ void getComponents(ComplexElement* e, Vector3 const& param, double_complex* comp int countNodes(ComplexElement* e) { - return countNodes(static_cast(e)); + return countNodes(e); } - +void getShapeValues(ComplexElement* e, Vector3 const& local, NewArray& values) +{ + getShapeValues(e,local,values); +} +void getShapeGrads(ComplexElement* e, Vector3 const& local, NewArray& grads) +{ + getShapeGrads(e,local,grads); +} +void getPackedNodes(ComplexElement* e, NewArray& values) +{ + FieldBase * f = e->getFieldBase(); + int cmps = f->countComponents(); + ComplexElementOf* element = static_cast*>(e); + element->getValues(values,cmps); +} } diff --git a/apf/apfComplexField.h b/apf/apfComplexField.h index bee7ed81f..2c22ca8c4 100644 --- a/apf/apfComplexField.h +++ b/apf/apfComplexField.h @@ -8,17 +8,6 @@ namespace apf { -class ComplexElement : public ElementT -{ -public: - ComplexElement(FieldBase* f, MeshEntity* e) - : ElementT(f,e) - { } - ComplexElement(FieldBase* f, VectorElement* p) - : ElementT(f,p) - { } -}; - class ComplexField : public FieldBase { public: @@ -30,6 +19,70 @@ class ComplexField : public FieldBase virtual void axpy(double_complex a, Field* x) = 0; }; +template +class ComplexFieldOf; + +template +void project(ComplexFieldOf* to, ComplexFieldOf* from) {} + +template +void axpy(double_complex a, ComplexFieldOf* x, ComplexFieldOf* y) {} + +template +class ComplexFieldOf::value>::type > : public ComplexField +{ +public: + virtual ~ComplexFieldOf() {} + void setNodeValue(MeshEntity* e, int node, T const value) + { + getData()->setNodeComponents( + e,node,reinterpret_cast(value)); + } + void getNodeValue(MeshEntity* e, int node, T* value) + { + getData()->getNodeComponents( + e,node,reinterpret_cast(value)); + } + void project(ComplexField* from) + { + apf::project(this,static_cast*>(from)); + } + void axpy(double_complex a, ComplexField* x) + { + apf::axpy(this,a,reinterpret_cast*>(x)); + } +}; + +template +class ComplexElementOf : public ComplexElement +{ +public: + ComplexElementOf(ComplexFieldOf* f, MeshEntity* e) + : Element(f,e) + { } + ComplexElementOf(ComplexFieldOf* f, VectorElement* p) + : Element(f,p) + { } + virtual ~ComplexElementOf() { } + T* getNodeValues() + { + return reinterpret_cast(&(this->nodeData[0])); + } + T getValue(Vector3 const& local) + { + T value[1]; + getComponents(local, reinterpret_cast(value)); + return value[0]; + } + void getValues(NewArray& 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: diff --git a/apf/apfComplexType.h b/apf/apfComplexType.h new file mode 100644 index 000000000..f3c437337 --- /dev/null +++ b/apf/apfComplexType.h @@ -0,0 +1,10 @@ +#ifdef C_COMPLEX + #include +#endif + +#define CXX_COMPLEX 1 +#ifdef CXX_COMPLEX + #include + using double_complex = std::complex; +#endif + diff --git a/apf/apfElement.cc b/apf/apfElement.cc index 193ca5f23..70e1e8c67 100644 --- a/apf/apfElement.cc +++ b/apf/apfElement.cc @@ -8,32 +8,39 @@ #include "apfElement.h" #include "apfShape.h" #include "apfMesh.h" +#include "apfComplexType.h" #include "apfVectorElement.h" namespace apf { -void ElementBase::init(FieldBase* f, MeshEntity* e, VectorElement* p) +template +ElementBase::ElementBase(FieldBase* f, MeshEntity* e) { - field = f; - mesh = f->getMesh(); - entity = e; - shape = f->getShape()->getEntityShape(mesh->getType(e)); - parent = p; - nen = shape->countNodes(); - nc = f->countComponents(); - getNodeData(); + init(f,e,0); } -ElementBase::ElementBase(FieldBase* f, MeshEntity* e) +template +ElementBase::ElementBase(FieldBase* f, VectorElement* p) { - init(f,e,0); + init(f,p->getEntity(),p); } -ElementBase::ElementBase(FieldBase* f, VectorElement* p) +template +void ElementBase::getGlobalGradients(Vector3 const& local, NewArray& globalGradients) { - init(f,p->getEntity(),p); + Matrix3x3 J; + parent->getJacobian(local,J); + Matrix3x3 jinv = getJacobianInverse(J, getDimension()); + NewArray 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; +template class ElementBase; + Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim) { switch (dim) { @@ -68,22 +75,4 @@ Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim) } } -int countNodes(ElementBase * e) -{ - return e->getShape()->countNodes(); -} - -void ElementBase::getGlobalGradients(Vector3 const& local, - NewArray& globalGradients) -{ - Matrix3x3 J; - parent->getJacobian(local,J); - Matrix3x3 jinv = getJacobianInverse(J, getDimension()); - NewArray localGradients; - shape->getLocalGradients(mesh, entity, local,localGradients); - globalGradients.allocate(nen); - for (int i=0; i < nen; ++i) - globalGradients[i] = jinv * localGradients[i]; -} - }//namespace apf diff --git a/apf/apfElement.h b/apf/apfElement.h index 797f91003..1dc17d424 100644 --- a/apf/apfElement.h +++ b/apf/apfElement.h @@ -15,9 +15,12 @@ namespace apf { +Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim); + class EntityShape; class VectorElement; +template class ElementBase { public: @@ -34,28 +37,6 @@ class ElementBase Mesh* getMesh() { return mesh; } EntityShape * getShape() { return shape; } FieldBase* getFieldBase() { return field; } -protected: - void init(FieldBase* f, MeshEntity* e, VectorElement* p); - virtual void getNodeData() = 0; - FieldBase* field; - Mesh* mesh; - MeshEntity* entity; - EntityShape* shape; - VectorElement* parent; - int nen; - int nc; -}; - -// not to be confused with ElementOf, which is about -// Scalar/Vector/Matrix node/value types -// this is about the underlying scalar type (real/complex) -template -class ElementT : public ElementBase -{ -public: - ElementT(FieldBase* f, MeshEntity* e) : ElementBase(f,e) { } - ElementT(FieldBase* f, VectorElement* p) : ElementBase(f,p) { } - virtual ~ElementT() {} void getComponents(Vector3 const& xi, T * c) { NewArray shapeValues; @@ -67,26 +48,49 @@ class ElementT : public ElementBase c[ci] += nodeData[ni * nc + ci] * shapeValues[ni]; } protected: - virtual void getNodeData() + 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*>(field->getData())->getElementData(entity,nodeData); } NewArray nodeData; + FieldBase* field; + Mesh* mesh; + MeshEntity* entity; + EntityShape* shape; + VectorElement* parent; + int nen; + int nc; }; -class Element : public ElementT + +template +int countNodes(ElementBase * e) { - public: - Element(FieldBase* f, MeshEntity* e) - : ElementT(f,e) - { } - Element(FieldBase* f, VectorElement* p) - : ElementT(f,p) - { } -}; + return e->getShape()->countNodes(); +} -Matrix3x3 getJacobianInverse(Matrix3x3 J, int dim); -int countNodes(ElementBase * e); +template +void getShapeValues(ElementBase * e, Vector3 const& local, NewArray& values) +{ + e->getShape()->getValues(e->getMesh(),e->getEntity(),local,values); +} + +template +void getShapeGrads(ElementBase * e, Vector3 const& local, NewArray& grads) +{ + e->getGlobalGradients(local,grads); +} }//namespace apf diff --git a/apf/apfElementOf.h b/apf/apfElementOf.h index 104826f94..923232331 100644 --- a/apf/apfElementOf.h +++ b/apf/apfElementOf.h @@ -14,8 +14,9 @@ namespace apf { -template ::value, bool> = 0 > +// std::enable_if_t::value, bool> = 0 > + +template class ElementOf : public Element { public: diff --git a/apf/apfField.h b/apf/apfField.h index 3659b44ce..6b5703871 100644 --- a/apf/apfField.h +++ b/apf/apfField.h @@ -10,10 +10,10 @@ #include #include "apfMesh.h" +#include "apfElementType.h" namespace apf { -class Element; class FieldData; class FieldBase diff --git a/apf/apfFieldOf.h b/apf/apfFieldOf.h index 94bad5d7b..765663f9e 100644 --- a/apf/apfFieldOf.h +++ b/apf/apfFieldOf.h @@ -13,10 +13,8 @@ namespace apf { -template -using enable = typename std::enable_if_t::value, bool>; -template +template class FieldOf; template @@ -26,8 +24,11 @@ void project(FieldOf* to, FieldOf* from); template void axpy(double a, FieldOf* x, FieldOf* y); +// template +// class FieldOf; +// typename std::enable_if::value>::type template -class FieldOf::value>::type > : public Field +class FieldOf : public Field { public: virtual ~FieldOf() {} @@ -43,11 +44,11 @@ class FieldOf::value>::typ } void project(Field* from) { - apf::project(this,static_cast*>(from)); + apf::project(this,static_cast*>(from)); } void axpy(double a, Field* x) { - apf::axpy(a,static_cast*>(x),this); + apf::axpy(a,static_cast*>(x),this); } }; diff --git a/apf/apfMesh.cc b/apf/apfMesh.cc index 912d92829..2a7dd0dbd 100644 --- a/apf/apfMesh.cc +++ b/apf/apfMesh.cc @@ -8,6 +8,7 @@ #include #include "apfCoordData.h" #include "apfVectorField.h" +#include "apfComplex.h" #include "apfComplexField.h" #include "apfShape.h" #include "apfNumbering.h" diff --git a/apf/apfMesh.h b/apf/apfMesh.h index 178b32932..ae9714bc8 100644 --- a/apf/apfMesh.h +++ b/apf/apfMesh.h @@ -14,7 +14,7 @@ #include #include #include -#include "apfComplex.h" +#include "apfComplexType.h" #include "apfVector.h" #include "apfDynamicArray.h" @@ -22,6 +22,7 @@ struct gmi_model; namespace apf { +class ComplexField; class FieldShape; class Field; class FieldBase; diff --git a/apf/apfZero.h b/apf/apfZero.h index 250ede35b..bfee2bfdc 100644 --- a/apf/apfZero.h +++ b/apf/apfZero.h @@ -2,6 +2,7 @@ #define APFZERO_H #include "apfField.h" +#include "apfNew.h" namespace apf { 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..bdd2e5e19 --- /dev/null +++ b/test/complex.cc @@ -0,0 +1,99 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +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 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] == double_complex(0,0) && r[1] == double_complex(0,0)); + 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 dofs1(cmps*nds1); + apf::NewArray 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/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 From b587a2590182acea5e1dfdbe31f91d6314fe05d6 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Tue, 9 Jul 2019 18:21:12 -0400 Subject: [PATCH 09/17] reverting accidental change of pField->apf::Field inside of pumi --- pumi/pumi_ghost.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pumi/pumi_ghost.cc b/pumi/pumi_ghost.cc index 5ada109ed..b4bb70091 100644 --- a/pumi/pumi_ghost.cc +++ b/pumi/pumi_ghost.cc @@ -400,7 +400,7 @@ void pumi_ghost_create(pMesh m, Ghosting* plan) std::vector frozen_fields; for (int i=0; icountFields(); ++i) { - apf::Field * f = m->getField(i); + pField * f = m->getField(i); if (isFrozen(f)) { frozen_fields.push_back(f); // turn field data from tag to array From 91f67a9221e5f1e4da3e5fe2e92add54275a8607 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Tue, 9 Jul 2019 18:23:12 -0400 Subject: [PATCH 10/17] adding apfComplexType.h to public interface --- apf/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/apf/CMakeLists.txt b/apf/CMakeLists.txt index 0c853935e..626d4cb4a 100644 --- a/apf/CMakeLists.txt +++ b/apf/CMakeLists.txt @@ -76,6 +76,7 @@ set(HEADERS apfFieldData.h apfNumberingClass.h apfComplex.h + apfComplexType.h ) # Add the apf library From 99423ccce83a9e925dba82bc72992e96a549d4f5 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Tue, 9 Jul 2019 18:25:17 -0400 Subject: [PATCH 11/17] adding new public headers to the pkg_tribits.cmake file --- apf/pkg_tribits.cmake | 2 ++ 1 file changed, 2 insertions(+) diff --git a/apf/pkg_tribits.cmake b/apf/pkg_tribits.cmake index a0c6cff0c..93a305d61 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 From 63744bea107b27b7abd0775126ddb1b682a49c09 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Tue, 9 Jul 2019 18:54:33 -0400 Subject: [PATCH 12/17] adding CMake option to determine compex type, C-complex is the default, had to change all occurences of the symbol 'I' in the codebase because the c-complex header defines a macro with that name, so we cannot use that symbol anywhere in the codebase --- CMakeLists.txt | 5 +++++ apf/CMakeLists.txt | 2 +- apf/apfComplexType.h | 10 ++++----- crv/crvQuality.cc | 50 ++++++++++++++++++++++---------------------- pumi/pumi_ghost.cc | 2 +- test/complex.cc | 3 ++- test/quality.cc | 4 ++-- 7 files changed, 40 insertions(+), 36 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index ae9496d6c..5dd5d58f4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -84,6 +84,11 @@ endif() set(Trilinos_PREFIX "" CACHE STRING "Trilinos installation directory") +option(CXX_COMPLEX "Use std::complex instead of double_complex" OFF) +if(CXX_COMPLEX) + add_definitions(-DCXX_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 626d4cb4a..9001e29e0 100644 --- a/apf/CMakeLists.txt +++ b/apf/CMakeLists.txt @@ -96,7 +96,7 @@ target_link_libraries(apf lion can mth - ) + ) scorec_export_library(apf) diff --git a/apf/apfComplexType.h b/apf/apfComplexType.h index f3c437337..fcfb0406e 100644 --- a/apf/apfComplexType.h +++ b/apf/apfComplexType.h @@ -1,9 +1,7 @@ -#ifdef C_COMPLEX - #include -#endif - -#define CXX_COMPLEX 1 -#ifdef CXX_COMPLEX +#ifndef CXX_COMPLEX +#include +typedef double _Complex double_complex; +#else #include using double_complex = std::complex; #endif 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& nodes, * */ static double Nijk(apf::NewArray& 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& 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& 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& nodes) static void getTriJacDetNodes(int P, apf::NewArray& elemNodes, apf::NewArray& 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& elemNodes, // apf::NewArray& nodes) @@ -581,11 +581,11 @@ double computeTriJacobianDetFromBezierFormulation(apf::Mesh* m, apf::NewArray 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/pumi/pumi_ghost.cc b/pumi/pumi_ghost.cc index b4bb70091..af59bfd63 100644 --- a/pumi/pumi_ghost.cc +++ b/pumi/pumi_ghost.cc @@ -400,7 +400,7 @@ void pumi_ghost_create(pMesh m, Ghosting* plan) std::vector frozen_fields; for (int i=0; icountFields(); ++i) { - pField * f = m->getField(i); + pField f = m->getField(i); if (isFrozen(f)) { frozen_fields.push_back(f); // turn field data from tag to array diff --git a/test/complex.cc b/test/complex.cc index bdd2e5e19..d4f9c6b07 100644 --- a/test/complex.cc +++ b/test/complex.cc @@ -21,6 +21,7 @@ int main(int ac, char** av) 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} }; { @@ -29,7 +30,7 @@ int main(int ac, char** av) while ((vert = m->iterate(it))) { apf::getComponents(f,vert,0,&r[0]); - PCU_ALWAYS_ASSERT(r[0] == double_complex(0,0) && r[1] == double_complex(0,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]); 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); From 0b151f87a2162e295e196214da1d600e6b5d0a0a Mon Sep 17 00:00:00 2001 From: William Tobin Date: Tue, 9 Jul 2019 20:12:31 -0400 Subject: [PATCH 13/17] forgot a file --- apf/apfElementType.h | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) create mode 100644 apf/apfElementType.h 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 ElementBase; +typedef ElementBase Element; +typedef ElementBase ComplexElement; + +}; + +#endif From e8c761b51e9462c04998aaba16df629be448f5b3 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Tue, 9 Jul 2019 20:34:17 -0400 Subject: [PATCH 14/17] removing reference to an unneeded file not in the repo --- apf/apfComplex.h | 1 - 1 file changed, 1 deletion(-) diff --git a/apf/apfComplex.h b/apf/apfComplex.h index 39cf9a06a..4b61497c3 100644 --- a/apf/apfComplex.h +++ b/apf/apfComplex.h @@ -2,7 +2,6 @@ #define APFCOMPLEX_H_ #include "apfComplexType.h" -#include "apfComplexElement.h" namespace apf { From 5624ba1dc2b2470dce1593a0bb83dc29c3d9fbb4 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Tue, 9 Jul 2019 21:07:17 -0400 Subject: [PATCH 15/17] fixing include issue caused by resolving unneeded intermediate file --- apf/apfComplex.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/apf/apfComplex.h b/apf/apfComplex.h index 4b61497c3..78301b253 100644 --- a/apf/apfComplex.h +++ b/apf/apfComplex.h @@ -2,11 +2,13 @@ #define APFCOMPLEX_H_ #include "apfComplexType.h" +#include "apfElementType.h" namespace apf { // forward decls for the interface +class ComplexField; class Mesh; class FieldShape; class MeshEntity; From 3aa1ebe14ccb83e876ca3eb3fc10ca8dc9bd628d Mon Sep 17 00:00:00 2001 From: William Tobin Date: Wed, 10 Jul 2019 09:16:40 -0400 Subject: [PATCH 16/17] addining a public header to the install list as it is needed downstream --- apf/CMakeLists.txt | 1 + apf/pkg_tribits.cmake | 1 + 2 files changed, 2 insertions(+) diff --git a/apf/CMakeLists.txt b/apf/CMakeLists.txt index 9001e29e0..83911cda2 100644 --- a/apf/CMakeLists.txt +++ b/apf/CMakeLists.txt @@ -77,6 +77,7 @@ set(HEADERS apfNumberingClass.h apfComplex.h apfComplexType.h + apfElementType.h ) # Add the apf library diff --git a/apf/pkg_tribits.cmake b/apf/pkg_tribits.cmake index 93a305d61..e7de49943 100644 --- a/apf/pkg_tribits.cmake +++ b/apf/pkg_tribits.cmake @@ -75,6 +75,7 @@ set(APF_HEADERS apfField.h apfFieldData.h apfNumberingClass.h + apfElementType.h ) set(APF_SOURCES From a9ebda59d55e08f54ae856b2a4fe5962f5816441 Mon Sep 17 00:00:00 2001 From: William Tobin Date: Wed, 10 Jul 2019 10:07:48 -0400 Subject: [PATCH 17/17] removing commented line of code, propogating the complex decision downstream without needing to carry along a configuration flag --- CMakeLists.txt | 8 +++++--- apf/CMakeLists.txt | 3 +++ apf/{apfComplexType.h => apfComplexType.h.in} | 4 +++- apf/apfElementOf.h | 2 -- 4 files changed, 11 insertions(+), 6 deletions(-) rename apf/{apfComplexType.h => apfComplexType.h.in} (74%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5dd5d58f4..f9ea9d1d6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -84,9 +84,11 @@ endif() set(Trilinos_PREFIX "" CACHE STRING "Trilinos installation directory") -option(CXX_COMPLEX "Use std::complex instead of double_complex" OFF) -if(CXX_COMPLEX) - add_definitions(-DCXX_COMPLEX) +option(PUMI_CXX_COMPLEX "Use std::complex 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) diff --git a/apf/CMakeLists.txt b/apf/CMakeLists.txt index 83911cda2..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 diff --git a/apf/apfComplexType.h b/apf/apfComplexType.h.in similarity index 74% rename from apf/apfComplexType.h rename to apf/apfComplexType.h.in index fcfb0406e..6ea989de1 100644 --- a/apf/apfComplexType.h +++ b/apf/apfComplexType.h.in @@ -1,4 +1,6 @@ -#ifndef CXX_COMPLEX +#define @PUMI_COMPLEX@ + +#ifndef PUMI_CXX_COMPLEX #include typedef double _Complex double_complex; #else diff --git a/apf/apfElementOf.h b/apf/apfElementOf.h index 923232331..fe0aabc75 100644 --- a/apf/apfElementOf.h +++ b/apf/apfElementOf.h @@ -14,8 +14,6 @@ namespace apf { -// std::enable_if_t::value, bool> = 0 > - template class ElementOf : public Element {