add function to update user field function

this is important if the apf::Function stores a field for computation.
It allows the user to update the apf::Function to use the proper field.

Author: jacobmerson

apf/apf.cc

@@ -463,6 +463,13 @@ Field* createUserField(Mesh* m, const char* name, int valueType, FieldShape* s,
   return makeField(m, name, valueType, 0, s, new UserData(f));
 }
 
+void updateUserField(Field* field, Function* newFunc)
+{
+  UserData* ud = dynamic_cast<UserData*>(field->getData());
+  // ud will be null if the data is not user data
+  if (ud) ud->setFunction(newFunc);
+}
+
 void copyData(Field* to, Field* from)
 {
   copyFieldData(from->getData(), to->getData());

apf/apf.h

@@ -668,14 +668,22 @@ struct Function
   virtual void eval(MeshEntity* e, double* result) = 0;
 };
 
-/** \brief Create a Field from a user's analytic function.
-  \details This field will use no memory, has values on all
-  nodes, and calls the user Function for all value queries.
-  Writing to this field does nothing.
-  */
+/* \brief Create a Field from a user's analytic function.
+ * \details This field will use no memory, has values on all
+ * nodes, and calls the user Function for all value queries.
+ * Writing to this field does nothing.
+ * \warning if you copy the mesh with apf::convert you may want to use
+ * apf::updateUserField to update function that this field refers to. This is
+ * extremely important if the analytic function you use references user fields.
+ */
 Field* createUserField(Mesh* m, const char* name, int valueType, FieldShape* s,
     Function* f);
 
+/* \brief update the analytic function from a user field
+ * \details this field updates the apf::Function which the UserField refers to
+ */
+void updateUserField(Field* field, Function* newFunc);
+
 /** \brief Compute a nodal gradient field from a nodal input field
   \details given a nodal field, compute approximate nodal gradient
   values by giving each node a volume-weighted average of the

apf/apfUserData.h

@@ -22,11 +22,14 @@ struct UserData : public FieldDataOf<double>
   void set(MeshEntity* e, double const* data);
   bool isFrozen();
   virtual FieldData* clone();
+  // using const * const gives an error on gcc/7.3.0 because the return is an
+  // r-value which cannot be modified anyways
+  Function const* getFunction() const { return function; }
+  void setFunction(Function* func) { function = func; }
+  private:
   Function* function;
 };
 
 }
 
 #endif
-
-

test/verify_convert.cc

@@ -1,60 +1,90 @@
 /*
  * this test verifies that the convert process properly clones the underlying
- * fields and numberings
+ * fields
  */
+#include <PCU.h>
 #include <apf.h>
+#include <apfConvert.h>
+#include <apfMDS.h>
+#include <apfMesh2.h>
 #include <gmi_null.h>
-#include <PCU.h>
 #include <pcu_util.h>
-#include <apfMesh2.h>
-#include <apfMDS.h>
-#include <apfConvert.h>
-
-apf::Mesh2* createEmptyMesh() {
+#include <cassert>
+apf::Mesh2* createEmptyMesh()
+{
   gmi_model* mdl = gmi_load(".null");
   return apf::makeEmptyMdsMesh(mdl, 3, false);
 }
-apf::Mesh2* createMesh() {
+apf::Mesh2* createMesh()
+{
   apf::Mesh2* m = createEmptyMesh();
   apf::MeshEntity* verts[2];
-  verts[0] = m->createVertex(NULL, apf::Vector3(0,0,0), apf::Vector3(0,0,0));
-  verts[1] = m->createVertex(NULL, apf::Vector3(1,0,0), apf::Vector3(1,0,0));
+  verts[0] =
+      m->createVertex(NULL, apf::Vector3(0, 0, 0), apf::Vector3(0, 0, 0));
+  verts[1] =
+      m->createVertex(NULL, apf::Vector3(1, 0, 0), apf::Vector3(1, 0, 0));
   m->createEntity(apf::Mesh::Type::EDGE, NULL, verts);
   return m;
 }
+class twox : public apf::Function {
+  private:
+  apf::Field* x;
 
-int main(int argc, char* argv[]) {
-  MPI_Init(&argc,&argv);
+  public:
+  virtual void eval(apf::MeshEntity* e, double* result)
+  {
+    result[0] = 2 * apf::getScalar(x, e, 0);
+  }
+  twox(apf::Field* x) : x(x) {}
+};
+int main(int argc, char* argv[])
+{
+  MPI_Init(&argc, &argv);
   PCU_Comm_Init();
   gmi_register_null();
   apf::Mesh* m1 = createMesh();
   apf::Mesh2* m2 = createEmptyMesh();
   // create field on m1
   apf::Field* f = apf::createLagrangeField(m1, "field1", apf::SCALAR, 1);
+  apf::Function* func = new twox(f);
+  apf::Field* uf =
+      apf::createUserField(m1, "ufield1", apf::SCALAR, apf::getShape(f), func);
   // loop over all vertices in mesh and set values
   int count = 1;
   apf::MeshIterator* it = m1->begin(0);
-  while(apf::MeshEntity * vert = m1->iterate(it)) {
+  while (apf::MeshEntity* vert = m1->iterate(it)) {
     apf::setScalar(f, vert, 0, count++);
   }
   m1->end(it);
+  // verify the user field works on mesh 1
+  it = m1->begin(0);
+  while (apf::MeshEntity* vert = m1->iterate(it)) {
+    double val = apf::getScalar(f, vert, 0);
+    double uval = apf::getScalar(uf, vert, 0);
+    if (!(std::abs(uval - 2 * double(val)) < 1E-15)) return 1;
+  }
+  m1->end(it);
   // copy m1 to m2
-  apf::convert(m1, m2); 
-  // read values from m2 and verify that they match values on m1
-  // here we assume that the nodes are ordered in the same way for
-  // both meshes
+  apf::convert(m1, m2);
   apf::Field* f2 = m2->findField("field1");
+  apf::Field* uf2 = m2->findField("ufield1");
+  // update the user field to reference the field in mesh 2
+  apf::updateUserField(uf2, new twox(f2));
   count = 1;
   it = m2->begin(0);
-  while(apf::MeshEntity * vert = m2->iterate(it) ) {
+  while (apf::MeshEntity* vert = m2->iterate(it)) {
     double val = apf::getScalar(f2, vert, 0);
-    if(!(std::abs(val-count) < 1E-15))
-      return 1;
+    double uval = apf::getScalar(uf2, vert, 0);
+    assert(std::abs(val - count) < 1E-15);
+    assert(std::abs(uval - 2 * double(count)) < 1E-15);
+    apf::setScalar(f2, vert, 0, 18.0);
+    // make sure that the function updated properly
+    uval = apf::getScalar(uf2, vert, 0);
+    assert(std::abs(uval - 36.0) < 1E-15);
     ++count;
   }
   m2->end(it);
-
-
+  delete func;
   m1->destroyNative();
   m2->destroyNative();
   apf::destroyMesh(m1);