MeshSmoothingOpFactory.hpp

Go to the documentation of this file.

#ifndef __MESH_SMOOTHING_OP_FACTORY_HPP__
#define __MESH_SMOOTHING_OP_FACTORY_HPP__
 
#ifndef WITH_ADOL_C
  #error "MoFEM need to be compiled with ADOL-C"
#endif
 
#include <MoFEM.hpp>
 
#include <SurfaceSlidingConstrains.hpp>
#include <MatOps.hpp>
#include <MatElastic.hpp>
 
#include <algorithm>
#include <limits>
#include <vector>
 
namespace MeshSmoothingOps {
 
struct Tags {
  int smoothingTag = 1;
  int surfaceTag = 2;
  int edgeTag = 3;
};
 
struct MinQualityOp;
 
struct Config {
  std::string geometryField = "MESH_NODE_POSITIONS";
  std::string meshReferenceField = "NONE";
  std::string lambdaSurfaceField = "LAMBDA_SURFACE";
  std::string lambdaEdgeField = "LAMBDA_EDGE";
  bool enableSurfaceSliding = true;
  bool enableEdgeSliding = false;
  bool fixEdgeBlocks = true;
  MatOps::VolumeLengthQualityType qualityType = MatOps::BARRIER_AND_QUALITY;
  double gamma = 0.0;
  Tags tags;
};
 
struct Handles {
  MoFEM::Interface *m_field_ptr;
  boost::shared_ptr<MatOps::PhysicalEquations> volQualityPhysicalPtr;
  boost::shared_ptr<ForcesAndSourcesCore> feSmootherRhs;
  boost::shared_ptr<ForcesAndSourcesCore> feSmootherLhs;
  boost::shared_ptr<SurfaceSlidingConstrains> surfaceConstraint;
  boost::shared_ptr<SurfaceSlidingConstrains::DriverElementOrientation>
      surfaceOrientation;
  boost::shared_ptr<EdgeSlidingConstrains> edgeConstraint;
  boost::shared_ptr<MoFEM::VolumeElementForcesAndSourcesCore> minQualityFe;
  boost::shared_ptr<double> minQualityValue;
};
 
struct MinQualityOp : MoFEM::ForcesAndSourcesCore::UserDataOperator {
 
  double *minQualityPtr;
 
  MinQualityOp(const std::string &field_name, double *min_quality_ptr)
      : MoFEM::ForcesAndSourcesCore::UserDataOperator(field_name,
                                                      UserDataOperator::OPROW),
        minQualityPtr(min_quality_ptr) {}
 
  MoFEMErrorCode doWork(int, EntityType type,
                        EntitiesFieldData::EntData &data) {
    MoFEMFunctionBeginHot;
    if (type != MBVERTEX)
      MoFEMFunctionReturnHot(0);
    const double q =
        MoFEM::Tools::volumeLengthQuality(&*data.getFieldData().begin());
    *minQualityPtr = std::min(*minQualityPtr, q);
    MoFEMFunctionReturnHot(0);
  }
};
 
inline MoFEMErrorCode createOperators(MoFEM::Interface &m_field,
                                      const Config &cfg,
                                      const Range &fixed_vertices,
                                      const Range &edge_entities,
                                      Handles &handles) {
  MoFEMFunctionBegin;
 
  handles.m_field_ptr = &m_field;
 
  using DomainEle = MoFEM::VolumeElementForcesAndSourcesCore;
  using DomainEleOp = DomainEle::UserDataOperator;
  handles.feSmootherRhs = boost::make_shared<DomainEle>(m_field);
  handles.feSmootherLhs = boost::make_shared<DomainEle>(m_field);
 
  handles.volQualityPhysicalPtr = MatOps::createMatOpsPhysicalEquationsPtr<
      MatOps::PIOLA_RESPONSE::VOLUMELENGTHQUALITY, MatOps::MODEL_3D>(
      MatOps::createMatOpsDataPtr(),
      MatOps::MatOpsTagsRegistry::setTagName("VolumeLength",
                                             cfg.tags.smoothingTag));
  CHKERR handles.volQualityPhysicalPtr->getOptions(&m_field);
  CHKERR handles.volQualityPhysicalPtr->recordTape();
 
  using DomainEleOp =
      MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator;
  auto &rhs_ops = handles.feSmootherRhs->getOpPtrVector();
  rhs_ops.clear();
  CHKERR MatOps::OpMaterialFactory<MatOps::PIOLA_RESPONSE, MatOps::MODEL_3D>::
      template opRhsFactory<PETSC, GAUSS, DomainEleOp>(
          m_field, rhs_ops, cfg.geometryField, handles.volQualityPhysicalPtr);
 
  auto &lhs_ops = handles.feSmootherLhs->getOpPtrVector();
  lhs_ops.clear();
  CHKERR MatOps::OpMaterialFactory<MatOps::PIOLA_RESPONSE, MatOps::MODEL_3D>::
      template opLhsFactory<PETSC, GAUSS, DomainEleOp>(
          m_field, lhs_ops, cfg.geometryField, handles.volQualityPhysicalPtr);
 
  handles.minQualityFe =
      boost::make_shared<MoFEM::VolumeElementForcesAndSourcesCore>(m_field);
  handles.minQualityValue =
      boost::make_shared<double>(std::numeric_limits<double>::max());
 
  handles.minQualityFe->getOpPtrVector().push_back(
      new MinQualityOp(cfg.geometryField, handles.minQualityValue.get()));
 
  if (!fixed_vertices.empty()) {
    auto *prb_mng = m_field.getInterface<ProblemsManager>();
    auto *simple = m_field.getInterface<Simple>();
 
    CHKERR prb_mng->removeDofsOnEntities(
        simple->getProblemName(), cfg.geometryField, fixed_vertices, 0, 3);
 
    if (cfg.enableSurfaceSliding &&
        m_field.check_field(cfg.lambdaSurfaceField)) {
      CHKERR prb_mng->removeDofsOnEntities(simple->getProblemName(),
                                           cfg.lambdaSurfaceField,
                                           fixed_vertices, 0, 1);
    }
 
    if (cfg.enableEdgeSliding && !cfg.fixEdgeBlocks &&
        m_field.check_field(cfg.lambdaEdgeField)) {
      CHKERR prb_mng->removeDofsOnEntities(
          simple->getProblemName(), cfg.lambdaEdgeField, fixed_vertices, 0, 1);
    }
  }
 
  if (cfg.enableSurfaceSliding) {
    handles.surfaceOrientation = boost::make_shared<
        SurfaceSlidingConstrains::DriverElementOrientation>();
    handles.surfaceConstraint = boost::shared_ptr<SurfaceSlidingConstrains>(
        handles.surfaceOrientation,
        new SurfaceSlidingConstrains(m_field, *handles.surfaceOrientation));
    CHKERR handles.surfaceConstraint->setOperators(
        cfg.tags.surfaceTag, cfg.lambdaSurfaceField, cfg.geometryField);
  }
 
  if (cfg.enableEdgeSliding && !edge_entities.empty() && !cfg.fixEdgeBlocks) {
    Range tets;
    CHKERR m_field.get_moab().get_entities_by_type(0, MBTET, tets);
    Range skin_faces;
    Skinner skin(&m_field.get_moab());
    CHKERR skin.find_skin(0, tets, false, skin_faces);
    if (skin_faces.empty()) {
      SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
              "No skin faces found for edge constrains");
    }
    handles.edgeConstraint = boost::make_shared<EdgeSlidingConstrains>(m_field);
    CHKERR handles.edgeConstraint->setOperators(cfg.tags.edgeTag, edge_entities,
                                                skin_faces, cfg.lambdaEdgeField,
                                                cfg.geometryField);
  }
 
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode addToSnesDM(DM dm, const Config &cfg,
                                  const Handles &handles) {
  MoFEMFunctionBegin;
  boost::shared_ptr<ForcesAndSourcesCore> null;
 
  CHKERR DMMoFEMSNESSetFunction(dm, "SMOOTHING", handles.feSmootherRhs, null,
                                null);
  CHKERR DMMoFEMSNESSetJacobian(dm, "SMOOTHING", handles.feSmootherLhs, null,
                                null);
 
  if (cfg.enableSurfaceSliding && handles.surfaceConstraint) {
    CHKERR DMMoFEMSNESSetFunction(
        dm, "SURFACE_SLIDING", handles.surfaceConstraint->feRhsPtr, null, null);
    CHKERR DMMoFEMSNESSetJacobian(
        dm, "SURFACE_SLIDING", handles.surfaceConstraint->feLhsPtr, null, null);
  }
 
  if (cfg.enableEdgeSliding && !cfg.fixEdgeBlocks && handles.edgeConstraint) {
    CHKERR DMMoFEMSNESSetFunction(dm, "EDGE_SLIDING",
                                  handles.edgeConstraint->feRhsPtr, null, null);
    CHKERR DMMoFEMSNESSetJacobian(dm, "EDGE_SLIDING",
                                  handles.edgeConstraint->feLhsPtr, null, null);
  }
 
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode evaluateMinQuality(DM dm,
                                         const std::string &domain_fe_name,
                                         const Handles &handles,
                                         double &min_quality) {
  MoFEMFunctionBegin;
 
  if (!handles.minQualityFe || !handles.minQualityValue) {
    SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
            "Min quality operators are not initialised");
  }
 
  *handles.minQualityValue = std::numeric_limits<double>::max();
  CHKERR DMoFEMLoopFiniteElements(dm, domain_fe_name.c_str(),
                                  handles.minQualityFe);
 
  MoFEM::Interface *m_field_ptr = nullptr;
  CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
  double local_min = *handles.minQualityValue;
  CHKERR MPI_Allreduce(&local_min, &min_quality, 1, MPI_DOUBLE, MPI_MIN,
                       m_field_ptr->get_comm());
 
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode updateBarrierGamma(Handles &handles, double gamma) {
  MoFEMFunctionBegin;
  char gamma_value[64];
  CHKERR PetscSNPrintf(gamma_value, sizeof(gamma_value), "%.17g", gamma);
  CHKERR PetscOptionsSetValue(NULL, "-volume_length_gamma", gamma_value);
  CHKERR handles.volQualityPhysicalPtr->getOptions(handles.m_field_ptr);
  CHKERR handles.volQualityPhysicalPtr->recordTape();
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode createSnes(DM dm, SmartPetscObj<SNES> &snes) {
  MoFEMFunctionBegin;
  MoFEM::Interface *m_field_ptr = nullptr;
  CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
  snes = createSNES(m_field_ptr->get_comm());
  CHKERR SNESSetFromOptions(snes);
  CHKERR SNESSetDM(snes, dm);
  CHKERR SNESSetUp(snes);
  MoFEMFunctionReturn(0);
}
 
inline MoFEMErrorCode
solveSnes(DM dm, SNES snes,
          const std::vector<std::string> &zero_vertex_fields = {}) {
  MoFEMFunctionBegin;
  if (!zero_vertex_fields.empty()) {
    MoFEM::Interface *m_field_ptr = nullptr;
    CHKERR DMoFEMGetInterfacePtr(dm, &m_field_ptr);
    auto field_blas = m_field_ptr->getInterface<FieldBlas>();
    for (const auto &field_name : zero_vertex_fields) {
      CHKERR field_blas->setField(0, MBVERTEX, field_name);
    }
  }
 
  auto f = createDMVector(dm);
  auto d = vectorDuplicate(f);
  CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR SNESSolve(snes, f, d);
  CHKERR DMoFEMMeshToGlobalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
  MoFEMFunctionReturn(0);
}
 
} // namespace MeshSmoothingOps
 
#endif // __MESH_SMOOTHING_OP_FACTORY_HPP__