Spr based error estimation

m3dc1_scorec/api/m3dc1_scorec.h

@@ -45,6 +45,7 @@ enum m3dc1_matrix_status { /*0*/ M3DC1_NOT_FIXED=0,
 bool m3dc1_double_isequal(double A, double B);
 
 int m3dc1_scorec_init();
+int m3dc1_scorec_verbosity(int*);
 int m3dc1_scorec_finalize();
 
 /** plane functions */
@@ -133,7 +134,8 @@ int m3dc1_field_getnewid (FieldID* /*out*/field_id);
 // is num_dofs input or output?
 // *value_type is either M3DC1_REAL or M3DC1_COMPLEX
 int m3dc1_field_create (FieldID* /*in*/ field_id, const char* /* in */ field_name, int* num_values, int* value_type, int* num_dofs_per_value);
-int m3dc1_field_delete (FieldID* /*in*/ field_id); 
+int m3dc1_mark_for_solutiontransfer (FieldID* /*in*/ field_id);
+int m3dc1_field_delete (FieldID* /*in*/ field_id);
 
 int m3dc1_field_getinfo(FieldID* /*in*/ field_id, char* /* out*/ field_name, int* num_values, int* value_type, int* total_num_dof);
 
@@ -215,6 +217,8 @@ int m3dc1_matrix_print(int* matrix_id);
 #endif // #ifdef M3DC1_PETSC
 
 // adaptation
+int m3dc1_spr_then_adapt (int * fieldId, int * index, int * ts,
+    double * ar, double * max_size, int * refine_level, int * coarsen_level, bool* update);
 int adapt_by_field (int * fieldId, double* psi0, double * psil);
 int set_adapt_p (double * pp);
 int adapt_by_error_field (double * errorField, double * errorAimed, int* max_node, int* option); // option 0: local error control; 1 global

m3dc1_scorec/api/name_convert.h

@@ -1,4 +1,5 @@
 #define m3dc1_scorec_init m3dc1_domain_init_
+#define m3dc1_scorec_verbosity m3dc1_domain_verbosity_
 #define m3dc1_scorec_finalize m3dc1_domain_finalize_
 #define m3dc1_plane_setnum m3dc1_plane_setnum_
 #define m3dc1_plane_getnum m3dc1_plane_getnum_
@@ -16,6 +17,7 @@
 #define m3dc1_mesh_load m3dc1_mesh_load_
 #define m3dc1_mesh_load_3d m3dc1_mesh_load_3d_
 #define m3dc1_mesh_build3d m3dc1_mesh_build3d_
+#define m3dc1_spr_then_adapt m3dc1_spr_then_adapt_
 #define m3dc1_mesh_adapt m3dc1_mesh_adapt_
 #define m3dc1_mesh_write m3dc1_mesh_write_
 #define m3dc1_mesh_getnument m3dc1_mesh_getnument_
@@ -59,6 +61,7 @@
 #define m3dc1_region_getoriginalface m3dc1_region_getoriginalface_
 #define m3dc1_field_getnewid m3dc1_field_genid_
 #define m3dc1_field_create m3dc1_field_create_
+#define m3dc1_mark_for_solutiontransfer m3dc1_mark_for_solutiontransfer_
 #define m3dc1_field_delete m3dc1_field_delete_
 #define m3dc1_field_exist m3dc1_field_exist_
 #define m3dc1_field_sync m3dc1_field_sync_

m3dc1_scorec/include/m3dc1_field.h

@@ -17,19 +17,23 @@ void group_complex_dof (apf::Field* field, int option);
 class m3dc1_field
 {
 public:
-  m3dc1_field (int i, apf::Field* f, int n, int t, int ndof): id(i), field(f), num_value(n), value_type(t),dof_per_value(ndof) {}
+  m3dc1_field (int i, apf::Field* f, int n, int t, int ndof): id(i), field(f), num_value(n), value_type(t),dof_per_value(ndof) {transfer = false;}
   ~m3dc1_field() {}
   apf::Field* get_field() { return field; }
+  void set_field(apf::Field* f) { field = f; }
   int get_id() { return id; }
   int get_num_value() { return num_value; }
   int get_value_type() { return value_type; }
   int get_dof_per_value() {return dof_per_value;}
+  bool should_transfer() {return transfer;}
+  void mark_for_solutiontransfer() {transfer = true;}
 private:
   int id;
   apf::Field* field; // name and #dofs are available from apf::Field
   int num_value;
   int value_type;
   int dof_per_value;
+  bool transfer;
 };
 
 void synchronize_field(apf::Field* f);

m3dc1_scorec/include/m3dc1_mesh.h

@@ -58,7 +58,9 @@ class m3dc1_mesh
   void restore3D();
 
   void build3d(int num_field, int* field_id, int* num_dofs_per_value);
-  void initialize(); 
+  void initialize();
+  void rebuildPointersOnNonMasterPlane(std::vector<std::vector<apf::Field*>>& pFields,
+              std::vector<apf::Field*>& zFields);
   void set_mcount(); // fill in # local, own, global mesh entity count
   void update_partbdry(apf::MeshEntity** remote_vertices, apf::MeshEntity** remote_edges, 
               apf::MeshEntity** remote_faces, std::vector<apf::MeshEntity*>& btw_plane_edges, 

m3dc1_scorec/src/m3dc1_adapt.cc

@@ -14,7 +14,9 @@
 #include <stdio.h>
 #include <stdlib.h>
 #include <PCU.h>
+#include <pcu_util.h>
 #include "apfMDS.h"
+#include "apfField.h"
 #include "ReducedQuinticImplicit.h"
 #include "m3dc1_slnTransfer.h"
 #include "apfShape.h" // getLagrange
@@ -115,7 +117,7 @@ void m3dc1_mesh::remove3D()
         mesh->destroy(e);
       mesh->end(ent_it);
     }
-  
+
   mesh->acceptChanges();
   changeMdsDimension(mesh, 2);
 
@@ -125,6 +127,187 @@ void m3dc1_mesh::remove3D()
     std::cout<<"\n*** Wedges and non-master 2D planes removed ***\n";
 }
 
+// this static function only used in m3dc1_mesh::rebuildPointersOnNonMasterPlane
+static int get_id_in_container(
+    std::map<FieldID, m3dc1_field*>& fcontainer,
+    apf::Field* f)
+{
+  typedef std::map<FieldID, m3dc1_field*> fct;
+
+  int id = -1;
+  for (fct::iterator it = fcontainer.begin(); it != fcontainer.end(); ++it)
+  {
+    m3dc1_field* mf = it->second;
+    if (f == mf->get_field())
+    {
+      id = it->first;
+      break;
+    }
+  }
+  return id;
+}
+// this will rebuild the mesh data-structure on non-master plane after remove3D
+//
+// Notes
+// (A) this is done to ensure "restore3D" behaves the same as "build3d".
+// More specifically, the order of entity creation and iteration remains the same
+// on all non-master planes.
+//
+// (B) the field pointers passed in the arrays will be updated for non-master planes
+void m3dc1_mesh::rebuildPointersOnNonMasterPlane(
+    std::vector<std::vector<apf::Field*>>& pFields, // multi-plane fields
+    std::vector<apf::Field*>& zFields)             // fields that need to be zero-d
+{
+  if (!(m3dc1_model::instance()->num_plane)) // if 2D do nothing
+    return;
+
+  // store the names so they can be updated later on
+  std::vector<std::vector<std::string>> pFieldNames; pFieldNames.clear();
+  std::vector<std::string>              zFieldNames; zFieldNames.clear();
+  std::vector<std::string>              tempNames;
+
+  for (int i = 0; i < (int)pFields.size(); i++)
+  {
+    tempNames.clear();
+    for (int j = 0; j < pFields[i].size(); j++)
+    {
+      tempNames.push_back(std::string(apf::getName(pFields[i][j])));
+    }
+    PCU_ALWAYS_ASSERT(pFields[i].size() == tempNames.size());
+    pFieldNames.push_back(tempNames);
+  }
+  PCU_ALWAYS_ASSERT(pFieldNames.size() == pFields.size());
+
+  for (int i = 0; i < (int)zFields.size(); i++)
+    zFieldNames.push_back(std::string(apf::getName(zFields[i])));
+  PCU_ALWAYS_ASSERT(zFieldNames.size() == zFields.size());
+
+
+  apf::Mesh2* mesh = m3dc1_mesh::instance()->mesh;
+  std::vector<int> fncs; fncs.clear();
+  std::vector<int> ids; ids.clear();  // if the field is in the m3dc1_mesh::instance()->field_container hold the id here
+  std::vector<std::string> fnames; fnames.clear();
+  std::vector<apf::FieldShape*> fshapes; fshapes.clear();
+
+  for (int i = 0; i < mesh->countFields(); i++) {
+    apf::Field* f = mesh->getField(i);
+    fncs.push_back(f->countComponents());
+    fnames.push_back(std::string(apf::getName(f)));
+    fshapes.push_back(apf::getShape(f));
+    int id = get_id_in_container(*m3dc1_mesh::instance()->field_container, f);
+    ids.push_back(id);
+  }
+
+  PCU_Barrier();
+
+  // now remove everything on non-master planes, and recreate them
+  // Removal Phase
+  // =============
+  // a-fields
+  // b-numberings
+  // c-tags
+  // d-internal mesh data-structure by calling destroyNative
+  // =============
+  // Recreate Phase
+  // a-create empty mds meshes on non-master planes
+  // b-set local_entid_tag, own_partid_tag, num_global_adj_node_tag, num_own_adj_node_tag to NULL
+  if (m3dc1_model::instance()->local_planeid != 0)
+  {
+    // manually delete all the fields/numberings and associated tags
+    while ( mesh->countFields() )
+    {
+      apf::Field* f = mesh->getField(0);
+      mesh->removeField(f);
+      apf::destroyField(f);
+    }
+
+    while ( mesh->countNumberings() )
+    {
+      apf::Numbering* n = mesh->getNumbering(0);
+      mesh->removeNumbering(n);
+      apf::destroyNumbering(n);
+    }
+
+    apf::DynamicArray<apf::MeshTag*> tags;
+    mesh->getTags(tags);
+    for (int i=0; i<tags.getSize(); i++)
+    {
+      if (mesh->findTag("norm_curv")==tags[i]) continue;
+      for (int idim=0; idim<4; idim++)
+        apf::removeTagFromDimension(mesh, tags[i], idim);
+      mesh->destroyTag(tags[i]);
+    }
+
+    // destroy the native mesh
+    apf::disownMdsModel(mesh);
+    mesh->destroyNative();
+    apf::destroyMesh(mesh);
+    m3dc1_mesh::instance()->mesh = pumi_mesh_create(pumi::instance()->model, 2, false);
+    mesh = m3dc1_mesh::instance()->mesh;
+
+    local_entid_tag = NULL;
+    own_partid_tag = NULL;
+    num_global_adj_node_tag = NULL;
+    num_own_adj_node_tag = NULL;
+  }
+
+  PCU_Barrier();
+
+
+  // update the field pointers in m3dc1_mesh::field_container
+  if (m3dc1_model::instance()->local_planeid != 0)
+  {
+    for (int i = 0; i < (int)fncs.size(); i++)
+      {
+	apf::Field* newf = apf::createPackedField(mesh, fnames[i].c_str(), fncs[i], fshapes[i]);
+	apf::zeroField(newf);
+	if (ids[i] > -1)
+	{
+	  m3dc1_field* mf = (*m3dc1_mesh::instance()->field_container)[ids[i]];
+	  mf->set_field(newf);
+	}
+      }
+  }
+
+  // update the field pointers in pFields
+  if (m3dc1_model::instance()->local_planeid != 0)
+  {
+    for (int i = 0; i < (int)pFieldNames.size(); i++)
+    {
+      for (int j = 0; j < pFieldNames[i].size(); j++)
+      {
+	apf::Field* f = mesh->findField(pFieldNames[i][j].c_str());
+	PCU_ALWAYS_ASSERT(f);
+	pFields[i][j] = f;
+      }
+    }
+  }
+
+
+  // update the field pointers in zFields
+  if (m3dc1_model::instance()->local_planeid != 0)
+  {
+    for (int i = 0; i < (int)zFieldNames.size(); i++)
+    {
+      apf::Field* f = mesh->findField(zFieldNames[i].c_str());
+      PCU_ALWAYS_ASSERT(f);
+      zFields[i] = f;
+    }
+  }
+
+
+  // since we have removed the Linear numbering for non-master planes, we add it beck here
+  apf::Numbering* linnumbering = mesh->findNumbering("Linear");
+  if (!linnumbering)
+    apf::createNumbering(mesh, "Linear", mesh->getShape(), 1);
+
+  mesh->acceptChanges();
+  set_mcount();
+
+  if (!PCU_Comm_Self())
+    std::cout<<"\n*** Data structures have been re-initialized on non-master planes ***\n";
+}
+
 void compute_size_and_frame_fields(apf::Mesh2* m, double* size_1, double* size_2, 
      double* angle, apf::Field* sizefield, apf::Field* framefield)
 {
@@ -356,12 +539,14 @@ void adapt_mesh (int field_id_h1, int field_id_h2, double* dir)
 
   ReducedQuinticImplicit shape;
   ReducedQuinticTransfer slnTransfer(mesh,fields, &shape);
+  ma::Input* in = ma::makeAdvanced(ma::configure(mesh, size_field, frame_field, &slnTransfer));
+/*
 #ifdef OLDMA
   ma::Input* in = ma::configure(mesh, size_field, frame_field, &slnTransfer);
 #else
   ma::Input* in = ma::makeAdvanced(ma::configure(mesh, size_field, frame_field, &slnTransfer));
 #endif
-
+*/	
   in->shouldSnap = 0;
   in->shouldTransferParametric = 0;
   in->shouldRunMidZoltan = 1;
@@ -1222,24 +1407,24 @@ void create_localid(apf::Mesh2* mesh, int dim)
 void m3dc1_mesh::restore3D()
 // *********************************************************
 {
-    // compute the fields to copy from master plane
+    /* // compute the fields to copy from master plane */
     int num_field = 0;
     int* field_id = NULL;
     int* num_dofs_per_value = NULL;
 
-    if (!m3dc1_model::instance()->local_planeid && m3dc1_mesh::instance()->field_container)
-    {
-      num_field = m3dc1_mesh::instance()->field_container->size();
-      cout<<__func__<<": #fields to copy from master to non-master plane "<<num_field<<"\n";
-      field_id = new int [num_field];
-      num_dofs_per_value = new int [num_field];
-      for (int i=0; i<num_field; ++i)
-      {
-        m3dc1_field * mf = (*(m3dc1_mesh::instance()->field_container))[*field_id];
-        field_id[i] = i;
-        num_dofs_per_value[i] = mf->get_dof_per_value();
-      }
-    }
+    /* if (!m3dc1_model::instance()->local_planeid && m3dc1_mesh::instance()->field_container) */
+    /* { */
+    /*   num_field = m3dc1_mesh::instance()->field_container->size(); */
+    /*   cout<<__func__<<": #fields to copy from master to non-master plane "<<num_field<<"\n"; */
+    /*   field_id = new int [num_field]; */
+    /*   num_dofs_per_value = new int [num_field]; */
+    /*   for (int i=0; i<num_field; ++i) */
+    /*   { */
+    /*     m3dc1_field * mf = (*(m3dc1_mesh::instance()->field_container))[*field_id]; */
+    /*     field_id[i] = i; */
+    /*     num_dofs_per_value[i] = mf->get_dof_per_value(); */
+    /*   } */
+    /* } */
 
   int local_partid=PCU_Comm_Self();
 
@@ -1487,9 +1672,9 @@ void m3dc1_mesh::restore3D()
   apf::Numbering* local_n = mesh->findNumbering(mesh->getShape()->getName());
   if (local_n) destroyNumbering(local_n);
 
-  // FIXME: re-create the field and copy field data on master process group to non-master
-  for (int i=0; i<num_field; ++i)
-    update_field(field_id[i], num_dofs_per_value[i], num_local_vtx, remote_vertices);
+  /* // FIXME: re-create the field and copy field data on master process group to non-master */
+  /* for (int i=0; i<num_field; ++i) */
+  /*   update_field(field_id[i], num_dofs_per_value[i], num_local_vtx, remote_vertices); */
 
   // clear temp memory
   delete [] remote_vertices;