Bezier Field Solution Transfer

crv/CMakeLists.txt

@@ -6,6 +6,8 @@ set(SOURCES
   crvBezier.cc
   crvBezierPoints.cc
   crvBezierShapes.cc
+  crvBezierFields.cc
+  crvBezierSolutionTransfer.cc
   crvBlended.cc
   crvCurveMesh.cc
   crvElevation.cc

crv/crv.h

@@ -11,6 +11,7 @@
 #include "apfMesh2.h"
 #include "apfShape.h"
 #include <ma.h>
+#include <maSolutionTransfer.h>
 #include <mth.h>
 #include <stdio.h>
 #include <vector>
@@ -22,6 +23,9 @@
   * \brief the curving functions are contained in this namespace */
 namespace crv {
 
+// forward declaration of the crv::Adapt
+class Adapt;
+
 /** \brief actually 1 greater than max order */
 static unsigned const MAX_ORDER = 19;
 
@@ -37,6 +41,13 @@ int getBlendingOrder(const int type);
 /** \brief count invalid elements of the mesh */
 int countNumberInvalidElements(apf::Mesh2* m);
 
+/** \ brief converts field f to Bezier entity wise */
+void convertInterpolationFieldPoints(apf::MeshEntity* e,
+    apf::Field* f, int n, int ne, apf::NewArray<double> &c);
+
+/** \brief converts field f, which is interpolating to Bezier */
+void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f);
+
 /** \brief Base Mesh curving object
   \details P is the order, S is the space dimension,
   different from the mesh dimension, used to distinguish between planar 2D
@@ -203,6 +214,10 @@ void writeInterpolationPointVtuFiles(apf::Mesh* m, const char* prefix);
 int getTriNodeIndex(int P, int i, int j);
 int getTetNodeIndex(int P, int i, int j, int k);
 
+/** \brief adds bezier solution transfers */
+ma::SolutionTransfer* setBezierSolutionTransfers(
+    const std::vector<apf::Field*>& fields, crv::Adapt* a);
+
 /** \brief crv fail function */
 void fail(const char* why) __attribute__((noreturn));
 

crv/crvAdapt.cc

@@ -194,6 +194,16 @@ static void flagCleaner(crv::Adapt* a)
   }
 }
 
+static void getAllBezierFields(ma::Mesh* m, std::vector<apf::Field*>& fields)
+{
+  for (int i = 0; i < m->countFields(); i++) {
+    apf::FieldShape* fs = apf::getShape(m->getField(i));
+    std::string name = fs->getName();
+    if (name == std::string("Bezier"))
+      fields.push_back(m->getField(i));
+  }
+}
+
 void adapt(ma::Input* in)
 {
   std::string name = in->mesh->getShape()->getName();
@@ -205,6 +215,14 @@ void adapt(ma::Input* in)
   double t0 = PCU_Time();
   ma::validateInput(in);
   Adapt* a = new Adapt(in);
+
+  // Setting up bezier field transfer for all fields with Bezier shapes
+  // This is not the cleanest way of doing this!
+  std::vector<apf::Field*> allFields;
+  getAllBezierFields(a->mesh, allFields);
+  in->solutionTransfer = crv::setBezierSolutionTransfers(allFields, a);
+  a->solutionTransfer = in->solutionTransfer;
+
   ma::preBalance(a);
 
   fixInvalidElements(a);

crv/crvBezierFields.cc

@@ -0,0 +1,136 @@
+/*
+ * Copyright 2015 Scientific Computation Research Center
+ *
+ * This work is open source software, licensed under the terms of the
+ * BSD license as described in the LICENSE file in the top-level directory.
+ */
+
+#include "crv.h"
+#include "crvBezier.h"
+#include "crvBezierShapes.h"
+#include "crvMath.h"
+#include "crvShape.h"
+#include "crvTables.h"
+#include "crvQuality.h"
+#include <apfField.h>
+#include <lionPrint.h>
+#include <cstdlib>
+#include <mth_def.h>
+#include <iostream>
+#include <pcu_util.h>
+
+namespace crv {
+
+static void convertVectorFieldInterpolationPoints(int n, int ne,
+    apf::NewArray<apf::Vector3>& nodes,
+    apf::NewArray<double>& c,
+    apf::NewArray<apf::Vector3>& newNodes){
+
+  for(int i = 0; i < ne; ++i)
+    newNodes[i].zero();
+
+  for( int i = 0; i < ne; ++i)
+    for( int j = 0; j < n; ++j)
+      newNodes[i] += nodes[j]*c[i*n+j];
+}
+
+static void convertScalarFieldInterpolationPoints(int n, int ne,
+    apf::NewArray<double>& nodes,
+    apf::NewArray<double>& c,
+    apf::NewArray<double>& newNodes)
+{
+
+  for(int i = 0; i < ne; ++i)
+    newNodes[i] = 0.;
+
+  for( int i = 0; i < ne; ++i)
+    for( int j = 0; j < n; ++j)
+      newNodes[i] += nodes[j]*c[i*n+j];
+}
+
+void convertInterpolationFieldPoints(apf::MeshEntity* e,
+    apf::Field* f,
+    int n, int ne, apf::NewArray<double>& c)
+{
+
+  apf::NewArray<apf::Vector3> l, b(ne);
+  apf::NewArray<double> ls, bs(ne);
+  apf::MeshElement* me =
+      apf::createMeshElement(f->getMesh(),e);
+  apf::Element* elem =
+      apf::createElement(f,me);
+
+  if (apf::getValueType(f) == apf::VECTOR) {
+    apf::getVectorNodes(elem,l);
+    convertVectorFieldInterpolationPoints(n, ne, l, c, b);
+    for(int i = 0; i < ne; ++i)
+      apf::setVector(f, e, i, b[i]);
+  }
+  else if (apf::getValueType(f) == apf::SCALAR) {
+    apf::getScalarNodes(elem, ls);
+    convertScalarFieldInterpolationPoints(n, ne, ls, c, bs);
+    for(int i = 0; i < ne; ++i)
+      apf::setScalar(f, e, i, bs[i]);
+  }
+  else
+    printf("Field type not implemented\n");
+  apf::destroyElement(elem);
+  apf::destroyMeshElement(me);
+}
+
+void convertInterpolatingFieldToBezier(apf::Mesh2* m_mesh, apf::Field* f)
+{
+  apf::FieldShape * fs = apf::getShape(f);
+  int order = fs->getOrder();
+
+  //apf::Field* fnew = createField(
+  //    m_mesh, "copy_field", apf::getValueType(f), crv::getBezier(order));
+  //transferFields(m_mesh, f, fnew);
+
+  int md = m_mesh->getDimension();
+  //
+  int blendingOrder = getBlendingOrder(apf::Mesh::simplexTypes[md]);
+  //blendingOrder = 0;
+  // go downward, and convert interpolating to control points
+  int startDim = md - (blendingOrder > 0);
+
+  for(int d = startDim; d >= 1; --d){
+    if(!fs->hasNodesIn(d)) continue;
+    int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes();
+    int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]);
+    apf::NewArray<double> c;
+    getBezierTransformationCoefficients(order,
+        apf::Mesh::simplexTypes[d],c);
+    apf::MeshEntity* e;
+    apf::MeshIterator* it = m_mesh->begin(d);
+    while ((e = m_mesh->iterate(it))){
+      if(m_mesh->isOwned(e))
+        convertInterpolationFieldPoints(e,f,n,ne,c);
+    }
+    m_mesh->end(it);
+  }
+
+  for( int d = 2; d <= md; ++d){
+  std::cout<<" blendingorder of dimension "<<d <<" "<<
+  getBlendingOrder(apf::Mesh::simplexTypes[d])<<std::endl;
+    if(!fs->hasNodesIn(d) ||
+        !getBlendingOrder(apf::Mesh::simplexTypes[d])) continue;
+    int n = fs->getEntityShape(apf::Mesh::simplexTypes[d])->countNodes();
+    int ne = fs->countNodesOn(apf::Mesh::simplexTypes[d]);
+    apf::NewArray<double> c;
+    getInternalBezierTransformationCoefficients(m_mesh,order,1,
+        apf::Mesh::simplexTypes[d],c);
+    apf::MeshEntity* e;
+    apf::MeshIterator* it = m_mesh->begin(d);
+    while ((e = m_mesh->iterate(it))){
+      if(!isBoundaryEntity(m_mesh,e) && m_mesh->isOwned(e))
+        convertInterpolationFieldPoints(e,f,n-ne,ne,c);
+    }
+    m_mesh->end(it);
+  }
+
+  apf::synchronize(f);
+
+}
+
+}