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 <Smoother.hpp>
#include <SurfaceSlidingConstrains.hpp>
#include <VolumeLengthQuality.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;
  VolumeLengthQualityType qualityType = BARRIER_AND_QUALITY;
  double alpha = 1.0;
  double gamma = 0.0;
  Tags tags;
};
 
struct Handles {
  boost::shared_ptr<Smoother> smoother;
  boost::shared_ptr<Smoother::MyVolumeFE> feSmootherRhs;
  boost::shared_ptr<Smoother::MyVolumeFE> feSmootherLhs;
  boost::shared_ptr<VolumeLengthQuality<double>> qualityDouble;
  boost::shared_ptr<VolumeLengthQuality<adouble>> qualityAdouble;
  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.smoother = boost::make_shared<Smoother>(m_field);
  handles.qualityAdouble = boost::make_shared<VolumeLengthQuality<adouble>>(
      cfg.qualityType, cfg.alpha, cfg.gamma, cfg.geometryField);
  handles.qualityDouble = boost::make_shared<VolumeLengthQuality<double>>(
      cfg.qualityType, cfg.alpha, cfg.gamma, cfg.geometryField);
 
  Range tets;
  CHKERR m_field.get_moab().get_entities_by_type(0, MBTET, tets);
  handles.smoother->setOfBlocks[0].tEts.merge(tets);
  handles.smoother->setOfBlocks[0].materialDoublePtr = handles.qualityDouble;
  handles.smoother->setOfBlocks[0].materialAdoublePtr = handles.qualityAdouble;
 
  handles.smoother->commonData.spatialPositions = cfg.geometryField;
  handles.smoother->commonData.meshPositions = cfg.meshReferenceField;
 
  handles.feSmootherRhs = handles.smoother->feRhsPtr;
  handles.feSmootherLhs = handles.smoother->feLhsPtr;
 
  auto &rhs_ops = handles.feSmootherRhs->getOpPtrVector();
  rhs_ops.clear();
  rhs_ops.push_back(new OpCalculateVectorFieldGradient<3, 3>(
      cfg.geometryField, handles.smoother->commonData.matGradPtr));
  rhs_ops.push_back(new Smoother::OpJacobianSmoother(
      cfg.geometryField, handles.smoother->setOfBlocks.at(0),
      handles.smoother->commonData, cfg.tags.smoothingTag, false));
  rhs_ops.push_back(new Smoother::OpInternalForcePiola(
      cfg.geometryField, handles.smoother->commonData.matFirstPiolaStressPtr));
 
  auto &lhs_ops = handles.feSmootherLhs->getOpPtrVector();
  lhs_ops.clear();
  lhs_ops.push_back(new OpCalculateVectorFieldGradient<3, 3>(
      cfg.geometryField, handles.smoother->commonData.matGradPtr));
  lhs_ops.push_back(new Smoother::OpJacobianSmoother(
      cfg.geometryField, handles.smoother->setOfBlocks.at(0),
      handles.smoother->commonData, cfg.tags.smoothingTag, true));
  lhs_ops.push_back(new Smoother::OpKPiola(
      cfg.geometryField, cfg.geometryField, handles.smoother->commonData.matTangentPtr));
 
  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 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;
  if (handles.qualityDouble)
    handles.qualityDouble->gAmma = gamma;
  if (handles.qualityAdouble)
    handles.qualityAdouble->gAmma = gamma;
  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__