Update to take into account last commits in Xfiles

CMakeLists.txt

@@ -34,13 +34,17 @@ find_package(xTool          REQUIRED)
 find_package(xInterfaceAOMDGeneral  REQUIRED)
 find_package(xDistMesh              REQUIRED )
 find_package(Trellis                REQUIRED )
+find_package(xMapping               REQUIRED)
+find_package(xQuadrature            REQUIRED )
 find_package(xSolverBase             REQUIRED )
 find_package(xInterfaceLapack        REQUIRED )
 find_package(xInterfaceSuperLu       REQUIRED )
 find_package(xInterfaceTaucs         REQUIRED )
+find_package(xInterfaceParMetis     REQUIRED )
 find_package(xFEM                   REQUIRED )
 find_package(xExt                   REQUIRED )
-find_package(xGeom                   REQUIRED )
+find_package(xGeomTools             REQUIRED)
+find_package(xDistanceNearest       REQUIRED)
 find_package(xTLS                   REQUIRED )
 find_package(xDoubleCut             REQUIRED )     
 find_package(xLegacySimpleCut       REQUIRED)
@@ -58,13 +62,16 @@ list(APPEND EXTERNAL_INCLUDES
   ${xInterfaceAOMDGeneral_INCLUDE_DIR}
   ${xDistMesh_INCLUDE_DIR}
   ${Trellis_INCLUDE_DIR}
+  ${xQuadrature_INCLUDE_DIR}
   ${xSolverBase_INCLUDE_DIR}
   ${xInterfaceSuperLu_INCLUDE_DIR}
   ${xInterfaceTaucs_INCLUDE_DIR}
   ${xInterfaceLapack_INCLUDE_DIR}
+  ${xInterfaceParMetis_INCLUDE_DIR}
   ${xFEM_INCLUDE_DIR}
   ${xExt_INCLUDE_DIR}
-  ${xGeom_INCLUDE_DIR}
+  ${xGeomTools_INCLUDE_DIR}
+  ${xDistanceNearest_INCLUDE_DIR}
   ${xTLS_INCLUDE_DIR}
   ${PARMETIS_INCLUDE_DIR}
   src
@@ -77,10 +84,13 @@ list(APPEND EXTERNAL_LIBRARIES
   ${xLegacySimpleCut_LIBRARIES}
   ${xDoubleCut_LIBRARIES}
   ${xExt_LIBRARIES}
-  ${xGeom_LIBRARIES}
+  ${xGeomTools_LIBRARIES}
+  ${xDistanceNearest_LIBRARIES}
   ${xFEM_LIBRARIES}
   ${xExport_LIBRARIES}
   ${Trellis_LIBRARIES}
+  ${xMapping_LIBRARIES} 
+  ${xQuadrature_LIBRARIES}
   ${xInterfaceAOMDGeneral_LIBRARIES}
   ${xSolverBase_LIBRARIES}
   ${xInterfaceLapack_LIBRARIES}

src/FormulationQSRemeshAniso.h

@@ -15,8 +15,8 @@ class TLSSolver;
 #include "xLevelSet.h"
 #include "xMesh.h"
 
-using xfem::xMesh;
 using AOMD::mEntity;
+using xfem::xMesh;
 
 class QSFormulRemeshAniso : public QSFormulation
 {
@@ -25,8 +25,8 @@ class QSFormulRemeshAniso : public QSFormulation
                        std::function<void(xMesh &)>);
    virtual ~QSFormulRemeshAniso();
 
-   //    void setBCDefinition(std::vector<std::tuple<std::function<bool (const Trellis_Util::mPoint&) >, int, int> >& BCDEF);
-   //    std::vector<std::tuple<std::function<bool (const Trellis_Util::mPoint&) >, int, int> > BCDEF;
+   //    void setBCDefinition(std::vector<std::tuple<std::function<bool (const xtensor::xPoint&) >, int, int> >& BCDEF);
+   //    std::vector<std::tuple<std::function<bool (const xtensor::xPoint&) >, int, int> > BCDEF;
    xfem::xLevelSet lsetAdapt;
 
    std::function<void(xMesh &)> tagger;
@@ -57,11 +57,11 @@ class xEvalOnOtherMesh : public xfem::xEval<T>
    void operator()(const xfem::xGeomElem *geo_appro, const xfem::xGeomElem *geo_integ, T &resu) const
    {
       // On recupere les coordonnees globales
-      Trellis_Util::mPoint xyz = geo_integ->getXYZ();  // On recupere les coordonnees du point considere
+      xtensor::xPoint xyz = geo_integ->getXYZ();  // On recupere les coordonnees du point considere
 
       // Recharge de l'element d'approx =============
-      std::set<AOMD::mEntity *> elts;     // On cree un "set" d'elements
-      meshInit.locateElement(xyz, elts);  // On cherche tous les elements qui contiennent xyz et on les met dans le set
+      std::vector<std::pair<const AOMD::mEntity *, const xtensor::xPoint>> elts =
+          meshInit.locateElement(xyz);  // On cherche tous les elements qui contiennent xyz et on les met dans le set
       // if(elts.begin()==elts.end()){std::cout<<"ERREUR dans le xEvalWiOnAll : Le noeud "<<xyz<<" n'est pas dans le
       // maillage\nContenu du-dit maillage :\n";meshInit.printAll();throw;}
 
@@ -74,8 +74,8 @@ class xEvalOnOtherMesh : public xfem::xEval<T>
       }
 
       // on extrait l'element parmis tous les elements
-      AOMD::mEntity *eltInit = 0;  // On cree un element unique
-      eltInit = *(elts.begin());   //...et on dit que c'est le premier element du set
+      AOMD::mEntity *eltInit = 0;                                      // On cree un element unique
+      eltInit = const_cast<AOMD::mEntity *>(((elts.begin())->first));  //...et on dit que c'est le premier element du set
 
 // Plus simple ici : les noeuds sont dupliques : xyz doit forcement etre place sur un noeud de meshInit
 // donc pas la peine de chercher la partition (qui n'existe pas, vue la maniere dont fonctionne le code ici)

src/FormulationQSSuperimposed.cc

-/* 
+/*
     This file is a part of eXlibris C++ Library
     under the GNU General Public License:
-    See the LICENSE.md files for terms and 
+    See the LICENSE.md files for terms and
     conditions.
 */
 
 #include "FormulationQSSuperimposed.h"
-#include "TLSGeom.h"
-#include "LinearSystem.h"
-
-#include "xAlgorithm.h"
-#include "xAssembler.h"
 
+#include "LinearSystem.h"
 #include "MaterialCommand.h"
 #include "NonUniformMaterialSensitivity.h"
+#include "TLSGeom.h"
+#include "xAlgorithm.h"
+#include "xAssembler.h"
 
 using namespace xfem;
 
-QSFormulationSuperimposed::QSFormulationSuperimposed(TLSGeom& geom, TLSSolver& tls_solver,
-                                                     const xData& data, const xParseData& parse_data,
-                                                     PreProcessing& pre_pro, PostProcessing& post_pro,
+QSFormulationSuperimposed::QSFormulationSuperimposed(TLSGeom& geom, TLSSolver& tls_solver, const xData& data,
+                                                     const xParseData& parse_data, PreProcessing& pre_pro,
+                                                     PostProcessing& post_pro,
                                                      const xEval<xtensor::xTensor2<>>& eval_strain_superimposed,
-                                                     const xEval<xtensor::xTensor2<>>& eval_stress_superimposed) :
-  QSFormulation(geom, tls_solver, data, parse_data, pre_pro, post_pro),
-  eval_strain_superimposed(eval_strain_superimposed),
-  eval_stress_superimposed(eval_stress_superimposed)
-{}
+                                                     const xEval<xtensor::xTensor2<>>& eval_stress_superimposed)
+    : QSFormulation(geom, tls_solver, data, parse_data, pre_pro, post_pro),
+      eval_strain_superimposed(eval_strain_superimposed),
+      eval_stress_superimposed(eval_stress_superimposed)
+{
+}
 
 QSFormulationSuperimposed::~QSFormulationSuperimposed() {}
 
-class MyNormal : public xEval<xtensor::xVector<>> {
-public:
-  void operator()(const xfem::xGeomElem* appro, const xfem::xGeomElem* integ, xtensor::xVector<>& vec) const {
-    auto pt=integ->getXYZ();
-    double theta=atan2(pt(1), pt(0));
-    vec(0)=cos(theta);
-    vec(1)=sin(theta);
-    vec(2)=0.;
-  }
+class MyNormal : public xEval<xtensor::xVector<>>
+{
+  public:
+   void operator()(const xfem::xGeomElem* appro, const xfem::xGeomElem* integ, xtensor::xVector<>& vec) const
+   {
+      auto pt = integ->getXYZ();
+      double theta = atan2(pt(1), pt(0));
+      vec(0) = cos(theta);
+      vec(1) = sin(theta);
+      vec(2) = 0.;
+   }
 };
 
-void QSFormulationSuperimposed::assembleNaturalEnv(LinearSystem& system) {
-  NonUniformMaterialSensitivity_c<xtensor::xTensor4<>> eval_stiffness("secant_matrix", variab_manager);
-  xEvalBinary<xtool::xMult<xtensor::xTensor4<>, xtensor::xTensor2<>, xtensor::xTensor2<>>> eval_stress(eval_stiffness, eval_strain_superimposed);
-  xFormLinearWithLoad<xGradOperator<xtensor::xSymmetrize<>>, xEval<xtensor::xTensor2<>> > form_linear(eval_stress);
-  xFilteredRegion<xIter, xEntityFilter> superimposed(geom.begin(), geom.end(), xAcceptOnZoneID(101));
-  system.setAssemblerCoeff(-1.);
-  Assemble(form_linear, system.getAssembler(), geom.getIntegRuleSmart(), disp_field, superimposed.begin(), superimposed.end());
-  system.setAssemblerCoeff(1.);
-  MyNormal eval_normal;
-  xEvalBinary<xtool::xMult<xtensor::xTensor2<>, xtensor::xVector<>, xtensor::xVector<>> > eval_traction(eval_stress_superimposed, eval_normal);
-  xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector<>> >, xEval<xtensor::xVector<>> > lin(eval_traction);
-  xClassRegion bc(&geom.getMesh(), 12, 1);
-  Assemble(lin, system.getAssembler(), geom.getIntegRuleBasic(bc_integ_order), disp_field, bc.begin(), bc.end(), xUpperAdjacency());
+void QSFormulationSuperimposed::assembleNaturalEnv(LinearSystem& system)
+{
+   NonUniformMaterialSensitivity_c<xtensor::xTensor4<>> eval_stiffness("secant_matrix", variab_manager);
+   xEvalBinary<xtool::xMult<xtensor::xTensor4<>, xtensor::xTensor2<>, xtensor::xTensor2<>>> eval_stress(eval_stiffness,
+                                                                                                        eval_strain_superimposed);
+   xFormLinearWithLoad<xGradOperator<xtensor::xSymmetrize<>>, xEval<xtensor::xTensor2<>>> form_linear(eval_stress);
+   xFilteredRegion<xIter, xEntityFilter> superimposed(geom.begin(), geom.end(), xAccept(101));
+   system.setAssemblerCoeff(-1.);
+   Assemble(form_linear, system.getAssembler(), geom.getIntegRuleSmart(), disp_field, superimposed.begin(), superimposed.end());
+   system.setAssemblerCoeff(1.);
+   MyNormal eval_normal;
+   xEvalBinary<xtool::xMult<xtensor::xTensor2<>, xtensor::xVector<>, xtensor::xVector<>>> eval_traction(eval_stress_superimposed,
+                                                                                                        eval_normal);
+   xFormLinearWithLoad<xValOperator<xtool::xIdentity<xtensor::xVector<>>>, xEval<xtensor::xVector<>>> lin(eval_traction);
+   xClassRegion bc(&geom.getMesh(), 12, 1);
+   Assemble(lin, system.getAssembler(), geom.getIntegRuleBasic(bc_integ_order), disp_field, bc.begin(), bc.end(),
+            xUpperAdjacency());
 }
 
-void QSFormulationSuperimposed::writeStrain(std::string name, std::string domain_name) {
-  OldAndCurrent_c::current();
-  xEvalGradField<xtensor::xSymmetrize<>> eval_strain_tilde(disp_field);
-  xEvalBinary<std::plus<xtensor::xTensor2<>> > eval_strain(eval_strain_tilde, eval_strain_superimposed);
-  SetMaterialVariablesVisitor_c<xEvalBinary<std::plus<xtensor::xTensor2<>>>> set_visitor(name, eval_strain);
-  UpdateMaterialVariablesVisitor_c update_visitor(name);
-  xeCompositeMaterialVariablesVisitor_c composite_visitor(set_visitor, update_visitor);
-  VisitMaterialVariablesCommand_c visit_command(composite_visitor, variab_manager);
-  ApplyCommandOnIntegrationRule(visit_command, geom.getIntegRuleSmart(), geom.begin(domain_name), geom.end(domain_name));
+void QSFormulationSuperimposed::writeStrain(std::string name, std::string domain_name)
+{
+   OldAndCurrent_c::current();
+   xEvalGradField<xtensor::xSymmetrize<>> eval_strain_tilde(disp_field);
+   xEvalBinary<std::plus<xtensor::xTensor2<>>> eval_strain(eval_strain_tilde, eval_strain_superimposed);
+   SetMaterialVariablesVisitor_c<xEvalBinary<std::plus<xtensor::xTensor2<>>>> set_visitor(name, eval_strain);
+   UpdateMaterialVariablesVisitor_c update_visitor(name);
+   xeCompositeMaterialVariablesVisitor_c composite_visitor(set_visitor, update_visitor);
+   VisitMaterialVariablesCommand_c visit_command(composite_visitor, variab_manager);
+   ApplyCommandOnIntegrationRule(visit_command, geom.getIntegRuleSmart(), geom.begin(domain_name), geom.end(domain_name));
 }

src/Import.cc

@@ -6,28 +6,26 @@
 */
 
 #include "Import.h"
+
+#include <fstream>
+
 #include "Export.h"
 #include "Formulation.h"
 #include "Observer.h"
 #include "TLSGeom.h"
 #include "TLSSolver.h"
 #include "Util.h"
-#include "options.h"
-
-#include "xAlgorithm.h"
-#include "xData.h"
-#include "xMesh.h"
-#include "xParseData.h"
-
-#include "workInProgress.h"
-
 #include "mAOMD.h"
 #include "mEntity.h"
 #include "mFace.h"
 #include "mTet.h"
 #include "mVertex.h"
-
-#include <fstream>
+#include "options.h"
+#include "workInProgress.h"
+#include "xAlgorithm.h"
+#include "xData.h"
+#include "xMesh.h"
+#include "xParseData.h"
 
 using namespace xfem;
 
@@ -219,9 +217,20 @@ int PreProcessing::getStep() const { return restart_step; }
 void PreProcessing::loadMesh(std::string mesh_name, xMesh& mesh)
 {
    std::string filename = getFilename(mesh_name);
-   AOMD::AOMD_Util::Instance()->import(filename.c_str(), &mesh);
-   mesh.modifyAllState();
-   mesh.bdryLinkSetup();
+   AOMD::AOMD_Util::Instance()->import(filename.c_str(), &(mesh.getMesh()));
+   mesh.getMesh().modifyState(3, 2, true);
+   mesh.getMesh().modifyState(3, 1, true);
+   mesh.getMesh().modifyState(3, 0, true);
+   mesh.getMesh().modifyState(2, 1, true);
+   mesh.getMesh().modifyState(2, 0, true);
+   mesh.getMesh().modifyState(1, 0, true);
+   mesh.getMesh().modifyState(0, 1, true);
+   mesh.getMesh().modifyState(0, 2, true);
+   mesh.getMesh().modifyState(0, 3, true);
+   mesh.getMesh().modifyState(1, 2, true);
+   mesh.getMesh().modifyState(1, 3, true);
+   mesh.getMesh().modifyState(2, 3, true);
+   mesh.getMesh().bdryLinkSetup();
    filename = "rm -f " + filename;
    system(filename.c_str());
 }
@@ -245,18 +254,18 @@ void PreProcessing::loadField(std::string field_name, xMesh& mesh, xField<>& fie
          double data[32];
          iss.read((char*)data, info[0] * sizeof(double));
 
-         AOMD::mVertex* n1 = mesh.getVertex(info[1]);
-         AOMD::mVertex* n2 = mesh.getVertex(info[2]);
-         AOMD::mVertex* n3 = mesh.getVertex(info[3]);
+         AOMD::mVertex* n1 = mesh.getMesh().getVertex(info[1]);
+         AOMD::mVertex* n2 = mesh.getMesh().getVertex(info[2]);
+         AOMD::mVertex* n3 = mesh.getMesh().getVertex(info[3]);
          AOMD::mEntity* e = nullptr;
          if (dim == 2)
          {
-            e = mesh.getTri(n1, n2, n3);
+            e = mesh.getMesh().getTri(n1, n2, n3);
          }
          else
          {
-            AOMD::mVertex* n4 = mesh.getVertex(info[4]);
-            e = mesh.getTet(n1, n2, n3, n4);
+            AOMD::mVertex* n4 = mesh.getMesh().getVertex(info[4]);
+            e = mesh.getMesh().getTet(n1, n2, n3, n4);
          }
          xFiniteElement FEM;
          FEM.setKeys(e, field.begin(), field.end());
@@ -297,18 +306,18 @@ void PreProcessing::loadDomain(std::string domain_name, xMesh& mesh, std::list<A
          int data[16];
          iss.read((char*)data, size * sizeof(int));
 
-         AOMD::mVertex* n1 = mesh.getVertex(data[0]);
-         AOMD::mVertex* n2 = mesh.getVertex(data[1]);
-         AOMD::mVertex* n3 = mesh.getVertex(data[2]);
+         AOMD::mVertex* n1 = mesh.getMesh().getVertex(data[0]);
+         AOMD::mVertex* n2 = mesh.getMesh().getVertex(data[1]);
+         AOMD::mVertex* n3 = mesh.getMesh().getVertex(data[2]);
          AOMD::mEntity* e = nullptr;
          if (dim == 2)
          {
-            e = mesh.getTri(n1, n2, n3);
+            e = mesh.getMesh().getTri(n1, n2, n3);
          }
          else
          {
-            AOMD::mVertex* n4 = mesh.getVertex(data[3]);
-            e = mesh.getTet(n1, n2, n3, n4);
+            AOMD::mVertex* n4 = mesh.getMesh().getVertex(data[3]);
+            e = mesh.getMesh().getTet(n1, n2, n3, n4);
          }
 
          list.push_back(e);

src/Integration.cc

-/* 
+/*
     This file is a part of eXlibris C++ Library
     under the GNU General Public License:
-    See the LICENSE.md files for terms and 
+    See the LICENSE.md files for terms and
     conditions.
 */
 
-#include "MeshGeneration.h"
 #include "Integration.h"
+
+#include "MeshGeneration.h"
+#include "mEntity.h"
 #include "xAttachableGP.h"
-#include "xIntegrationRuleStored.h"
 #include "xCSRVector.h"
 #include "xDenseMatrix.h"
+#include "xIntegrationRuleStored.h"
 #include "xLapackInterface.h"
-#include "mEntity.h"
 
 using namespace AOMD;
 using namespace xfem;
 
-IntegrationPointCreator::IntegrationPointCreator(mEntity* e_model,
-                                                 int degree,
-                                                 const unsigned int partition_tag,
-                                                 const unsigned int gp_tag) :
-  degree(degree),
-  rhs_degree(degree),
-  gauss_degree(1),
-  partition_tag(partition_tag),
-  gp_tag(gp_tag) {
-  assert(degree<5);
+IntegrationPointCreator::IntegrationPointCreator(mEntity* e_model, int degree, const unsigned int partition_tag,
+                                                 const unsigned int gp_tag)
+    : degree(degree), rhs_degree(degree), gauss_degree(1), partition_tag(partition_tag), gp_tag(gp_tag)
+{
+   assert(degree < 5);
 
-  const int nb_terms=(degree+1)*(degree+2)/2;
+   const int nb_terms = (degree + 1) * (degree + 2) / 2;
 
-  Trellis_Util::GaussIntegrator gauss_integrator(e_model);
-  int nb_gauss_points=gauss_integrator.nbIntegrationPoints(gauss_degree);
-  while(nb_gauss_points<nb_terms) {
-    ++gauss_degree;
-    nb_gauss_points=gauss_integrator.nbIntegrationPoints(gauss_degree);
-  }
-  if(nb_gauss_points>nb_terms) {
-    ++rhs_degree;
-  }
+   xquadrature::xGaussIntegrator gauss_integrator(e_model);
+   int nb_gauss_points = gauss_integrator.nbIntegrationPoints(gauss_degree);
+   while (nb_gauss_points < nb_terms)
+   {
+      ++gauss_degree;
+      nb_gauss_points = gauss_integrator.nbIntegrationPoints(gauss_degree);
+   }
+   if (nb_gauss_points > nb_terms)
+   {
+      ++rhs_degree;
+   }
 }
 
-void IntegrationPointCreator::operator()(mEntity* e_appro) const {
-  xfem::xGeomElem geo_appro(e_appro);
-  geo_appro.SetIntegrationPointNumberForDegree(gauss_degree);
-  int nb_gauss_points=geo_appro.GetNbIntegrationPoints();
-
-  std::vector<Trellis_Util::mPoint> points;
-  points.reserve(nb_gauss_points);
-  xlinalg::xDenseMatrix mat(nb_gauss_points, nb_gauss_points);
-  for(int k=0; k<nb_gauss_points; ++k) {
-    geo_appro.setUVW(k);
-    const Trellis_Util::mPoint& uvw=geo_appro.getUVW();
-    mat(0,k)=1.;
-    mat(1,k)=uvw(1);
-    mat(2,k)=uvw(0);
-    if(nb_gauss_points>3) {
-      mat(3,k)=uvw(1)*uvw(1);
-      mat(4,k)=uvw(0)*uvw(1);
-      mat(5,k)=uvw(0)*uvw(0);
-    }
-    if(nb_gauss_points>6) {
-      mat(6,k)=uvw(1)*uvw(1)*uvw(1);
-      mat(7,k)=uvw(0)*uvw(1)*uvw(1);
-      mat(8,k)=uvw(0)*uvw(0)*uvw(1);
-      mat(9,k)=uvw(0)*uvw(0)*uvw(0);
-      mat(10,k)=uvw(1)*uvw(1)*uvw(1)*uvw(1);
-      mat(11,k)=uvw(0)*uvw(1)*uvw(1)*uvw(1);
-    }
-    if(nb_gauss_points>12) {
-      mat(12,k)=uvw(0)*uvw(0)*uvw(1)*uvw(1);
-      mat(13,k)=uvw(0)*uvw(0)*uvw(0)*uvw(1);
-      mat(14,k)=uvw(0)*uvw(0)*uvw(0)*uvw(0);
-      mat(15,k)=uvw(1)*uvw(1)*uvw(1)*uvw(1)*uvw(1);
-    }
-    points.push_back(uvw);
-  }
+void IntegrationPointCreator::operator()(mEntity* e_appro) const
+{
+   xfem::xGeomElem geo_appro(e_appro);
+   geo_appro.SetIntegrationPointNumberForDegree(gauss_degree);
+   int nb_gauss_points = geo_appro.GetNbIntegrationPoints();
 
-  std::vector<mEntity*> partition;
-  AddFromPartition(partition_tag, partition)(e_appro);
-  xlinalg::xCSRVector vec(nb_gauss_points);
-  for(auto ep: partition) {
-    xfem::xGeomElem geo_integ(ep);
-    geo_integ.SetIntegrationPointNumberForDegree(rhs_degree);
-    const int nb=geo_integ.GetNbIntegrationPoints();
-    for(int k=0; k<nb; ++k) {
-      geo_integ.setUVW(k);
-      geo_appro.setUVWForXYZ(geo_integ.getXYZ());
-      const Trellis_Util::mPoint& uvw=geo_appro.getUVW();
-      const double coeff=geo_integ.GetWeight()*geo_integ.GetDetJac()/geo_appro.GetDetJac();
-      vec(0)+=coeff;
-      vec(1)+=uvw(1)*coeff;
-      vec(2)+=uvw(0)*coeff;
-      if(nb_gauss_points>3) {
-        vec(3)+=uvw(1)*uvw(1)*coeff;
-        vec(4)+=uvw(0)*uvw(1)*coeff;
-        vec(5)+=uvw(0)*uvw(0)*coeff;
+   std::vector<xtensor::xPoint> points;
+   points.reserve(nb_gauss_points);
+   xlinalg::xDenseMatrix mat(nb_gauss_points, nb_gauss_points);
+   for (int k = 0; k < nb_gauss_points; ++k)
+   {
+      geo_appro.setUVW(k);
+      const xtensor::xPoint& uvw = geo_appro.getUVW();
+      mat(0, k) = 1.;
+      mat(1, k) = uvw(1);
+      mat(2, k) = uvw(0);
+      if (nb_gauss_points > 3)
+      {
+         mat(3, k) = uvw(1) * uvw(1);
+         mat(4, k) = uvw(0) * uvw(1);
+         mat(5, k) = uvw(0) * uvw(0);
       }
-      if(nb_gauss_points>6) {
-        vec(6)+=uvw(1)*uvw(1)*uvw(1)*coeff;
-        vec(7)+=uvw(0)*uvw(1)*uvw(1)*coeff;
-        vec(8)+=uvw(0)*uvw(0)*uvw(1)*coeff;
-        vec(9)+=uvw(0)*uvw(0)*uvw(0)*coeff;
-        vec(10)+=uvw(1)*uvw(1)*uvw(1)*uvw(1)*coeff;
-        vec(11)+=uvw(0)*uvw(1)*uvw(1)*uvw(1)*coeff;
+      if (nb_gauss_points > 6)
+      {
+         mat(6, k) = uvw(1) * uvw(1) * uvw(1);
+         mat(7, k) = uvw(0) * uvw(1) * uvw(1);
+         mat(8, k) = uvw(0) * uvw(0) * uvw(1);
+         mat(9, k) = uvw(0) * uvw(0) * uvw(0);
+         mat(10, k) = uvw(1) * uvw(1) * uvw(1) * uvw(1);
+         mat(11, k) = uvw(0) * uvw(1) * uvw(1) * uvw(1);
       }
-      if(nb_gauss_points>12) {
-        vec(12)+=uvw(0)*uvw(0)*uvw(1)*uvw(1)*coeff;
-        vec(13)+=uvw(0)*uvw(0)*uvw(0)*uvw(1)*coeff;
-        vec(14)+=uvw(0)*uvw(0)*uvw(0)*uvw(0)*coeff;
-        vec(15)+=uvw(1)*uvw(1)*uvw(1)*uvw(1)*uvw(1)*coeff;
+      if (nb_gauss_points > 12)
+      {
+         mat(12, k) = uvw(0) * uvw(0) * uvw(1) * uvw(1);
+         mat(13, k) = uvw(0) * uvw(0) * uvw(0) * uvw(1);
+         mat(14, k) = uvw(0) * uvw(0) * uvw(0) * uvw(0);
+         mat(15, k) = uvw(1) * uvw(1) * uvw(1) * uvw(1) * uvw(1);
+      }
+      points.push_back(uvw);
+   }
+
+   std::vector<mEntity*> partition;
+   AddFromPartition(partition_tag, partition)(e_appro);
+   xlinalg::xCSRVector vec(nb_gauss_points);
+   for (auto ep : partition)
+   {
+      xfem::xGeomElem geo_integ(ep);
+      geo_integ.SetIntegrationPointNumberForDegree(rhs_degree);
+      const int nb = geo_integ.GetNbIntegrationPoints();
+      for (int k = 0; k < nb; ++k)
+      {
+         geo_integ.setUVW(k);
+         geo_appro.setUVWForXYZ(geo_integ.getXYZ());
+         const xtensor::xPoint& uvw = geo_appro.getUVW();
+         const double coeff = geo_integ.GetWeight() * geo_integ.GetDetJac() / geo_appro.GetDetJac();
+         vec(0) += coeff;
+         vec(1) += uvw(1) * coeff;
+         vec(2) += uvw(0) * coeff;
+         if (nb_gauss_points > 3)
+         {
+            vec(3) += uvw(1) * uvw(1) * coeff;
+            vec(4) += uvw(0) * uvw(1) * coeff;
+            vec(5) += uvw(0) * uvw(0) * coeff;
+         }
+         if (nb_gauss_points > 6)
+         {
+            vec(6) += uvw(1) * uvw(1) * uvw(1) * coeff;
+            vec(7) += uvw(0) * uvw(1) * uvw(1) * coeff;
+            vec(8) += uvw(0) * uvw(0) * uvw(1) * coeff;
+            vec(9) += uvw(0) * uvw(0) * uvw(0) * coeff;
+            vec(10) += uvw(1) * uvw(1) * uvw(1) * uvw(1) * coeff;
+            vec(11) += uvw(0) * uvw(1) * uvw(1) * uvw(1) * coeff;
+         }
+         if (nb_gauss_points > 12)
+         {
+            vec(12) += uvw(0) * uvw(0) * uvw(1) * uvw(1) * coeff;
+            vec(13) += uvw(0) * uvw(0) * uvw(0) * uvw(1) * coeff;
+            vec(14) += uvw(0) * uvw(0) * uvw(0) * uvw(0) * coeff;
+            vec(15) += uvw(1) * uvw(1) * uvw(1) * uvw(1) * uvw(1) * coeff;
+         }
       }
-    }
-  }
+   }
 
-  xlinalg::xCSRVector weights(nb_gauss_points);
-  solve(mat, vec, weights);
+   xlinalg::xCSRVector weights(nb_gauss_points);
+   solve(mat, vec, weights);
 
-  xAttachableGaussPoints* att=new xAttachableGaussPoints;
-  std::transform(points.begin(), points.end(),
-                 weights.begin(),
-                 std::back_inserter(att->gauss_points),
-                 std::make_pair<Trellis_Util::mPoint&, double&>);
-  e_appro->attachData(gp_tag, att);
+   xAttachableGaussPoints* att = new xAttachableGaussPoints;
+   std::transform(points.begin(), points.end(), weights.begin(), std::back_inserter(att->gauss_points),
+                  std::make_pair<xtensor::xPoint&, double&>);
+   e_appro->attachData(gp_tag, att);
 }

src/TLSGeom.cc

@@ -6,20 +6,19 @@
 */
 
 #include "TLSGeom.h"
+
 #include "Export.h"
 #include "Import.h"
 #include "MeshGeneration.h"
 #include "Observer.h"
-
 #include "xData.h"
 #include "xIntegrationRuleStored.h"
 #include "xLevelSetOperators.h"
 #include "xSubMesh.h"
 
 // xfem::distmesh
-#include "xLoadBalanceTools.h"
-
 #include "workInProgress.h"
+#include "xLoadBalanceTools.h"
 
 using namespace xfem;
 using AOMD::mEntity;
@@ -153,14 +152,14 @@ class DelUpdateAdj
    DelUpdateAdj(xMesh& mesh) : mesh(mesh) {}
    void operator()(mEntity* e)
    {
-      mesh.mMesh::del(e);
+      mesh.getMesh().mMesh::del(e);
       DelAdjacencies(e);
    }
 
   private:
    xMesh& mesh;
 };
-void ComputeBoundingBox(const xSubMesh& dom, Trellis_Util::mPoint& min, Trellis_Util::mPoint& max)
+void ComputeBoundingBox(const xSubMesh& dom, xtensor::xPoint& min, xtensor::xPoint& max)
 {
    auto it = dom.begin(0);
    auto ite = dom.end(0);
@@ -233,14 +232,14 @@ TLSGeom::TLSGeom(xData& data, const xParseData& parse_data, PreProcessing& pre_p
       // must be called only with one proc for now (tar split or not .... see import.cc)
       assert(xtool::workInProgress());
       // data.ReadMesh();
-      // comp_mesh=data.mesh;// TODO stored mesh as node position perturbed w.r.t. original mesh
+      // comp_mesh=data.getMesh();// TODO stored mesh as node position perturbed w.r.t. original mesh
       comp_mesh = new xMesh(world);
       pre_pro.loadMesh("comp", *comp_mesh);
    }
    else
    {
       data.ReadMesh(world);
-      comp_mesh = data.mesh;
+      comp_mesh = data.getMesh();
    }
 
    int nb_proc;
@@ -272,7 +271,7 @@ TLSGeom::TLSGeom(xData& data, const xParseData& parse_data, PreProcessing& pre_p
    xRegion comp(comp_mesh);
 
    h = xCFL::getLength(comp);
-   Trellis_Util::mPoint p_min, p_max;
+   xtensor::xPoint p_min, p_max;
    comp_mesh->compute_bounding_box(p_min, p_max);
    body_characteristic_length = std::max(p_max(0) - p_min(0), std::max(p_max(1) - p_min(1), p_max(2) - p_min(2)));
    std::cout << "characteristic element length " << h << std::endl;
@@ -479,7 +478,7 @@ void TLSGeom::buildNonlocal(const std::list<mEntity*>& delocalized_list)
 void TLSGeom::computeNonlocalBoundingBox() const
 {
    auto& sm_nonlocal = getMesh().getSubMesh("nonlocal");
-   Trellis_Util::mPoint min, max;
+   xtensor::xPoint min, max;
    ComputeBoundingBox(sm_nonlocal, min, max);
    int step = Observer::tell<int>("step");
    post_pro.exportOnTime("nonlocal_bb_min_x", step, min(0));
@@ -632,7 +631,7 @@ void TLSGeom::transferBndIntInNonlocal()
    auto& sm_int_in_nonlocal = mesh.getSubMesh("int_in_nonlocal");
 
    xAttached duplicated_parent(was_created_by_tag);
-   std::vector<mEntity *> todump_bnd, todump_int;
+   std::vector<mEntity*> todump_bnd, todump_int;
    todump_bnd.reserve(sm_int_duplicated.size(0));
    todump_int.reserve(sm_int_duplicated.size(0));
    for (auto it = sm_int_duplicated.begin(0); it != sm_int_duplicated.end(0); ++it)
@@ -1012,16 +1011,30 @@ void TLSGeom::buildCompMesh()
                  NodeDispatcher(comp_mesh, was_created_by_tag, is_crack_tag, partition_tag, duplicated_node_tag));
    std::for_each(sm_int_fds.begin(0), sm_int_fds.end(0),
                  NodeDuplicator(comp_mesh, was_created_by_tag, is_crack_tag, partition_tag, duplicated_node_tag));
-   std::for_each(sm_fds.begin(), sm_fds.end(), FaceDuplicator(comp_mesh, sm_duplicated, was_created_by_tag, is_crack_tag,
-                                                              partition_tag, duplicated_tag, duplicated_node_tag));
+   std::for_each(sm_fds.begin(), sm_fds.end(),
+                 FaceDuplicator(comp_mesh, sm_duplicated, was_created_by_tag, is_crack_tag, partition_tag, duplicated_tag,
+                                duplicated_node_tag));
    std::for_each(sm_fds.begin(), sm_fds.end(), DelUpdateAdj(comp_mesh));
    std::for_each(sm_int_fds.begin(1), sm_int_fds.end(1), DelUpdateAdj(comp_mesh));
    std::for_each(sm_int_fds.begin(0), sm_int_fds.end(0), DelUpdateAdj(comp_mesh));
 
-   comp_mesh.modifyAllState();  // TODO maybe improvable
-   comp_mesh.bdryLinkSetup();   // TODO maybe improvable
+   comp_mesh.getMesh().modifyState(3, 2, true);
+   comp_mesh.getMesh().modifyState(3, 1, true);
+   comp_mesh.getMesh().modifyState(3, 0, true);
+   comp_mesh.getMesh().modifyState(2, 1, true);
+   comp_mesh.getMesh().modifyState(2, 0, true);
+   comp_mesh.getMesh().modifyState(1, 0, true);
+   comp_mesh.getMesh().modifyState(0, 1, true);
+   comp_mesh.getMesh().modifyState(0, 2, true);
+   comp_mesh.getMesh().modifyState(0, 3, true);
+   comp_mesh.getMesh().modifyState(1, 2, true);
+   comp_mesh.getMesh().modifyState(1, 3, true);
+   comp_mesh.getMesh().modifyState(2, 3, true);
+
+   comp_mesh.getMesh().bdryLinkSetup();  // TODO maybe improvable
 
    sm_duplicated.modifyAllState();
+
    auto& sm_int_duplicated = comp_mesh.createSubMesh("int_duplicated");
    GetComplement(sm_duplicated, sm_bnd_fds, sm_int_duplicated, 0);
 

test/animesh_dam_notch/main.cc

-/* 
+/*
     This file is a part of eXlibris C++ Library
     under the GNU General Public License:
-    See the LICENSE.md files for terms and 
+    See the LICENSE.md files for terms and