Example Struct Reference

[Example] More...

Collaboration diagram for Example:

Classes
struct	BoundaryOp
struct	CommonData
	[Example] More...
struct	DynamicFirstOrderConsConstantTimeScale
struct	DynamicFirstOrderConsSinusTimeScale
struct	OpCalcSurfaceAverageTemperature
struct	OpError
struct	OpError< 1 >
struct	OpFirst
struct	OpFluxRhs
struct	OpRadiationLhs
struct	OpRadiationRhs
struct	OpRhs
struct	OpSecond
	[Operator] More...
struct	OpZero
	[Common data] More...
struct	ScaledTimeScale

Public Types
enum	{ VOL , COUNT }

Public Member Functions
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	[Run problem]
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Main driver function for the optimization process.

Static Public Attributes
static std::array< double, 2 >	meshVolumeAndCount = {0, 0}

Private Types
enum	BoundingBox {   CENTER_X = 0 , CENTER_Y , MAX_X , MAX_Y ,   MIN_X , MIN_Y , LAST_BB }

Private Member Functions
MoFEMErrorCode	setupProblem ()
	[Run problem]
MoFEMErrorCode	createCommonData ()
	[Set up problem]
MoFEMErrorCode	bC ()
	[Create common data]
MoFEMErrorCode	OPs ()
	[Boundary condition]
MoFEMErrorCode	tsSolve ()
MoFEMErrorCode	testOperators ()
	[Solve]
MoFEMErrorCode	readMesh ()
	[Run problem]
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	boundaryCondition ()
	[Set up problem]
MoFEMErrorCode	assembleSystem ()
	[Push operators to pipeline]
MoFEMErrorCode	solveSystem ()
	[Solve]
MoFEMErrorCode	outputResults ()
	[Solve]
MoFEMErrorCode	checkResults ()
	[Postprocess results]
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	setUp ()
	[Run all]
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	setFieldValues ()
	[Create common data]
MoFEMErrorCode	pushOperators ()
	[Set density distribution]
MoFEMErrorCode	integrateElements ()
	[Push operators to pipeline]
MoFEMErrorCode	postProcess ()
	[Integrate]
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
	[Set up problem]
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	assembleSystemIntensity ()
	[Calculate flux on boundary]
MoFEMErrorCode	assembleSystemFlux ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	calculateFlux (double &calc_flux)
	[Set up problem]
MoFEMErrorCode	outputResults (const int i)
	[Solve]
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	bC ()
MoFEMErrorCode	OPs ()
MoFEMErrorCode	kspSolve ()
	[Push operators to pipeline]
MoFEMErrorCode	postProcess ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
	Read mesh from file and setup meshsets.
MoFEMErrorCode	setupProblem ()
	Setup fields, approximation spaces and DOFs.
MoFEMErrorCode	setupAdJoint ()
	Setup adjoint fields and finite elements
MoFEMErrorCode	boundaryCondition ()
	Apply essential boundary conditions.
MoFEMErrorCode	topologyModes ()
	Compute topology optimization modes.
MoFEMErrorCode	assembleSystem ()
	Setup operators in finite element pipeline.
MoFEMErrorCode	solveElastic ()
	Solve forward elastic problem.
MoFEMErrorCode	postprocessElastic (int iter, SmartPetscObj< Vec > adjoint_vector=nullptr)
	Post-process and output results.
MoFEMErrorCode	calculateGradient (PetscReal *objective_function_value, Vec objective_function_gradient, Vec adjoint_vector)
	Calculate objective function gradient using adjoint method.

Static Private Member Functions
static std::pair< int, int >	getCoordsInImage (double x, double y)
static double	rhsSource (const double x, const double y, const double)
static double	lhsFlux (const double x, const double y, const double)
static int	integrationRule (int, int, int p_data)

Private Attributes
MoFEM::Interface &	mField
	Reference to MoFEM interface.
boost::shared_ptr< DomainEle >	reactionFe
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >	uXScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >	uYScatter
std::tuple< SmartPetscObj< Vec >, SmartPetscObj< VecScatter > >	uZScatter
Simple *	simpleInterface
boost::shared_ptr< std::vector< unsigned char > >	boundaryMarker
boost::shared_ptr< CommonData >	commonDataPtr
FieldApproximationBase	base
	Choice of finite element basis functions
FieldSpace	space
boost::shared_ptr< VectorDouble >	approxVals
boost::shared_ptr< MatrixDouble >	approxGradVals
Range	pinchNodes
boost::shared_ptr< MatrixDouble >	matDPtr
SmartPetscObj< Mat >	M
SmartPetscObj< Mat >	K
SmartPetscObj< EPS >	ePS
std::array< SmartPetscObj< Vec >, 6 >	rigidBodyMotion
boost::shared_ptr< FEMethod >	domianLhsFEPtr
boost::shared_ptr< FEMethod >	domianRhsFEPtr
boost::shared_ptr< MatrixDouble >	vectorFieldPtr = nullptr
	Field values at evaluation points.
int	fieldOrder = 2
	Polynomial order for approximation.
SmartPetscObj< KSP >	kspElastic
	Linear solver for elastic problem.
SmartPetscObj< DM >	adjointDM
	Data manager for adjoint problem.
boost::shared_ptr< ObjectiveFunctionData >	pythonPtr
	Interface to Python objective function.
std::vector< SmartPetscObj< Vec > >	modeVecs
	Topology mode vectors (design variables)
std::vector< std::array< double, 3 > >	modeCentroids
	Centroids of optimization blocks
std::vector< std::array< double, 6 > >	modeBBoxes
	Bounding boxes of optimization blocks.
SmartPetscObj< Vec >	initialGeometry
	Initial geometry field.

Static Private Attributes
static std::vector< double >	rZ
static std::vector< MatrixInt >	iI
static std::array< double, LAST_BB >	aveMaxMin
static int	focalIndex
static int	savitzkyGolayNormalisation
static const int *	savitzkyGolayWeights
static ApproxFieldFunction< FIELD_DIM >	approxFunction

Friends
struct	TSPrePostProc

Detailed Description

[Example]

Main class for topology optimization using adjoint method.

This class implements a complete topology optimization workflow:

1. Setup finite element problems for elasticity and adjoint analysis
2. Define topology modes for design parametrization
   
3. Use TAO optimization library to minimize Python-defined objective function
4. Compute sensitivities using adjoint method for efficient gradient calculation

The optimization process alternates between:

• Forward analysis: solve elastic problem for current geometry
• Adjoint analysis: solve adjoint problem for objective function gradients
• Geometry update: modify design using optimization algorithm (L-BFGS)

Examples

mofem/tutorials/adv-0/plastic.cpp, mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-1/integration.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/mix-1/phase.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-8/radiation.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4/shallow_wave.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 216 of file plastic.cpp.

Member Enumeration Documentation

anonymous enum

anonymous enum

Enumerator
VOL
COUNT

Definition at line 222 of file plastic.cpp.

{ VOL, COUNT };

BoundingBox

enum Example::BoundingBox

private

Enumerator
CENTER_X
CENTER_Y
MAX_X
MAX_Y
MIN_X
MIN_Y
LAST_BB

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 98 of file phase.cpp.

                   {
    CENTER_X = 0,
    CENTER_Y,
    MAX_X,
    MAX_Y,
    MIN_X,
    MIN_Y,
    LAST_BB
  };

Constructor & Destructor Documentation

Example() [1/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 218 of file plastic.cpp.

: mField(m_field) {}

Example() [2/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 420 of file dynamic_first_order_con_law.cpp.

: mField(m_field) {}

Example() [3/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 35 of file helmholtz.cpp.

: mField(m_field) {}

Example() [4/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 25 of file integration.cpp.

: mField(m_field) {}

Example() [5/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 50 of file plot_base.cpp.

: mField(m_field) {}

Example() [6/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 81 of file phase.cpp.

: mField(m_field) {}

Example() [7/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 53 of file approximaton.cpp.

: mField(m_field) {}

Example() [8/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 45 of file radiation.cpp.

: mField(m_field) {}

Example() [9/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 47 of file heat_method.cpp.

: mField(m_field) {}

Example() [10/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 49 of file eigen_elastic.cpp.

: mField(m_field) {}

Example() [11/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 52 of file nonlinear_dynamic_elastic.cpp.

: mField(m_field) {}

Example() [12/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 339 of file shallow_wave.cpp.

: mField(m_field) {}

Example() [13/13]

Example::Example ( MoFEM::Interface &  m_field )

inline

Definition at line 179 of file adjoint.cpp.

: mField(m_field) {}

Member Function Documentation

assembleSystem() [1/9]

MoFEMErrorCode Example::assembleSystem ( )

private

[Push operators to pipeline]

[Adjoint modes]

[Boundary condition]

[Applying essential BC]

[Push operators to pipeline]

[Integration rule]

[Integration rule]

[Push domain stiffness matrix to pipeline]

[Push domain stiffness matrix to pipeline]

[Push Body forces]

[Push Body forces]

[Push natural boundary conditions]

[Push natural boundary conditions]

Examples

mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4/shallow_wave.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 640 of file dynamic_first_order_con_law.cpp.

                                       {
  MoFEMFunctionBegin;
  auto get_body_force = [this](const double, const double, const double) {
    FTensor::Index<'i', SPACE_DIM> i;
    FTensor::Tensor1<double, SPACE_DIM> t_source;
    t_source(i) = 0.;
    t_source(0) = 0.1;
    t_source(1) = 1.;
    return t_source;
  };
 
  // specific time scaling
  auto get_time_scale = [this](const double time) {
    return sin(time * omega * M_PI);
  };
 
  auto apply_rhs = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
 
    // Calculate Gradient of velocity
    auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "V", mat_v_grad_ptr));
 
    auto gravity_vector_ptr = boost::make_shared<MatrixDouble>();
    gravity_vector_ptr->resize(SPACE_DIM, 1);
    auto set_body_force = [&]() {
      FTensor::Index<'i', SPACE_DIM> i;
      MoFEMFunctionBegin;
      auto t_force = getFTensor1FromMat<SPACE_DIM, 0>(*gravity_vector_ptr);
      double unit_weight = 0.;
      CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-unit_weight", &unit_weight,
                                 PETSC_NULLPTR);
      t_force(i) = 0;
      if (SPACE_DIM == 2) {
        t_force(1) = -unit_weight;
      } else if (SPACE_DIM == 3) {
        t_force(2) = unit_weight;
      }
      MoFEMFunctionReturn(0);
    };
 
    CHKERR set_body_force();
    pip.push_back(new OpBodyForce("V", gravity_vector_ptr,
                                  [](double, double, double) { return 1.; }));
 
    // Calculate unknown F
    auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(
        "F", mat_H_tensor_ptr));
 
    //  // Calculate F
    double tau = 0.2;
    CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-tau", &tau, PETSC_NULLPTR);
 
    double xi = 0.;
    CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-xi", &xi, PETSC_NULLPTR);
 
    // Calculate P stab
    auto one = [&](const double, const double, const double) {
      return 3. * bulk_modulus_K;
    };
    auto minus_one = [](const double, const double, const double) {
      return -1.;
    };
 
    auto mat_dot_F_tensor_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(
        "F_dot", mat_dot_F_tensor_ptr));
 
    // Calculate Gradient of Spatial Positions
    auto mat_x_grad_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "x_2", mat_x_grad_ptr));
 
    auto mat_F_tensor_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateDeformationGradient<SPACE_DIM, SPACE_DIM>(
        mat_F_tensor_ptr, mat_H_tensor_ptr));
 
    auto mat_F_stab_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateFStab<SPACE_DIM, SPACE_DIM>(
        mat_F_tensor_ptr, mat_F_stab_ptr, mat_dot_F_tensor_ptr, tau, xi,
        mat_x_grad_ptr, mat_v_grad_ptr));
 
    PetscBool is_linear_elasticity = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_linear_elasticity",
                               &is_linear_elasticity, PETSC_NULLPTR);
 
    auto mat_P_stab_ptr = boost::make_shared<MatrixDouble>();
    if (is_linear_elasticity) {
      pip.push_back(new OpCalculatePiola<SPACE_DIM, SPACE_DIM>(
          shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_stab_ptr,
          mat_F_stab_ptr));
    } else {
      auto inv_F = boost::make_shared<MatrixDouble>();
      auto det_ptr = boost::make_shared<VectorDouble>();
 
      pip.push_back(
          new OpInvertMatrix<SPACE_DIM>(mat_F_stab_ptr, det_ptr, inv_F));
 
      // OpCalculatePiolaIncompressibleNH
      pip.push_back(new OpCalculatePiolaIncompressibleNH<SPACE_DIM, SPACE_DIM>(
          shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_stab_ptr,
          mat_F_stab_ptr, inv_F, det_ptr));
    }
 
    pip.push_back(new OpGradTimesTensor2("V", mat_P_stab_ptr, minus_one));
    pip.push_back(new OpRhsTestPiola("F", mat_v_grad_ptr, one));
 
    MoFEMFunctionReturn(0);
  };
 
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  CHKERR apply_rhs(pipeline_mng->getOpDomainExplicitRhsPipeline());
 
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order;
  };
  CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);
 
  MoFEMFunctionReturn(0);
}

assembleSystem() [2/9]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [3/9]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [4/9]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [5/9]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [6/9]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [7/9]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [8/9]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [9/9]

MoFEMErrorCode Example::assembleSystem ( )

private

Setup operators in finite element pipeline.

assembleSystemFlux()

MoFEMErrorCode Example::assembleSystemFlux ( )

private

Examples

mofem/tutorials/mix-1/phase.cpp.

assembleSystemIntensity()

MoFEMErrorCode Example::assembleSystemIntensity ( )

private

[Calculate flux on boundary]

[Push operators to pipeline]

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 433 of file phase.cpp.

                                                {
  MoFEMFunctionBegin;
 
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  pipeline_mng->getDomainLhsFE().reset();
  pipeline_mng->getDomainRhsFE().reset();
  pipeline_mng->getBoundaryRhsFE().reset();
 
  auto rule_vol = [](int, int, int order) { return 2 * (order + 1); };
  pipeline_mng->setDomainLhsIntegrationRule(rule_vol);
  pipeline_mng->setDomainRhsIntegrationRule(rule_vol);
 
  auto det_ptr = boost::make_shared<VectorDouble>();
  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpCalculateHOJacForFace(jac_ptr));
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpMakeHdivFromHcurl());
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpSetInvJacHcurlFace(inv_jac_ptr));
  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetHOWeightsOnFace());
 
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpHdivHdiv("S", "S", lhsFlux));
  auto unity = []() { return 1; };
  pipeline_mng->getOpDomainLhsPipeline().push_back(
      new OpHdivU("S", "PHI", unity, true));
 
  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSetHOWeightsOnFace());
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpDomainSource("PHI", rhsSource));
 
  MoFEMFunctionReturn(0);
}

bC() [1/2]

MoFEMErrorCode Example::bC ( )

private

[Create common data]

[Boundary condition]

Examples

mofem/tutorials/adv-0/plastic.cpp, and mofem/tutorials/scl-8/radiation.cpp.

Definition at line 606 of file plastic.cpp.

                           {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_X",
                                           "U", 0, 0);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Y",
                                           "U", 1, 1);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), "REMOVE_Z",
                                           "U", 2, 2);
  CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(),
                                           "REMOVE_ALL", "U", 0, 3);
 
#ifdef ADD_CONTACT
  for (auto b : {"FIX_X", "REMOVE_X"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                             "SIGMA", 0, 0, false, true);
  for (auto b : {"FIX_Y", "REMOVE_Y"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                             "SIGMA", 1, 1, false, true);
  for (auto b : {"FIX_Z", "REMOVE_Z"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                             "SIGMA", 2, 2, false, true);
  for (auto b : {"FIX_ALL", "REMOVE_ALL"})
    CHKERR bc_mng->removeBlockDOFsOnEntities(simple->getProblemName(), b,
                                             "SIGMA", 0, 3, false, true);
  CHKERR bc_mng->removeBlockDOFsOnEntities(
      simple->getProblemName(), "NO_CONTACT", "SIGMA", 0, 3, false, true);
#endif
 
  CHKERR bc_mng->pushMarkDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "U");
 
  auto &bc_map = bc_mng->getBcMapByBlockName();
  for (auto bc : bc_map)
    MOFEM_LOG("PLASTICITY", Sev::verbose) << "Marker " << bc.first;
 
  MoFEMFunctionReturn(0);
}

bC() [2/2]

MoFEMErrorCode Example::bC ( )

private

boundaryCondition() [1/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

[Set up problem]

[Create common data]

[Boundary condition]

[Applying essential BC]

Examples

mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4/shallow_wave.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 535 of file dynamic_first_order_con_law.cpp.

                                          {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
  auto time_scale = boost::make_shared<TimeScale>();
 
  PetscBool sin_time_function = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-sin_time_function",
                             &sin_time_function, PETSC_NULLPTR);
 
  if (sin_time_function)
    time_scale = boost::make_shared<DynamicFirstOrderConsSinusTimeScale>();
  else
    time_scale = boost::make_shared<DynamicFirstOrderConsConstantTimeScale>();
 
  pipeline_mng->getBoundaryExplicitRhsFE().reset();
  CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
      pipeline_mng->getOpBoundaryExplicitRhsPipeline(), {NOSPACE}, "GEOMETRY");
 
  CHKERR BoundaryNaturalBC::AddFluxToPipeline<OpForce>::add(
      pipeline_mng->getOpBoundaryExplicitRhsPipeline(), mField, "V",
      {time_scale}, "FORCE", Sev::inform);
 
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order;
  };
 
  CHKERR pipeline_mng->setBoundaryExplicitRhsIntegrationRule(integration_rule);
  CHKERR pipeline_mng->setDomainExplicitRhsIntegrationRule(integration_rule);
 
  CHKERR bc_mng->removeBlockDOFsOnEntities<DisplacementCubitBcData>(
      simple->getProblemName(), "V");
 
  auto get_pre_proc_hook = [&]() {
    return EssentialPreProc<DisplacementCubitBcData>(
        mField, pipeline_mng->getDomainExplicitRhsFE(), {time_scale});
  };
  pipeline_mng->getDomainExplicitRhsFE()->preProcessHook = get_pre_proc_hook();
 
  MoFEMFunctionReturn(0);
}

boundaryCondition() [2/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [3/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [4/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [5/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [6/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [7/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [8/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [9/9]

MoFEMErrorCode Example::boundaryCondition ( )

private

Apply essential boundary conditions.

calculateFlux()

MoFEMErrorCode Example::calculateFlux ( double &  calc_flux )

private

[Set up problem]

[Calculate flux on boundary]

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 402 of file phase.cpp.

                                                       {
  MoFEMFunctionBegin;
  auto pipeline_mng = mField.getInterface<PipelineManager>();
 
  pipeline_mng->getDomainLhsFE().reset();
  pipeline_mng->getDomainRhsFE().reset();
  pipeline_mng->getBoundaryRhsFE().reset();
 
  auto rule_vol = [](int, int, int order) { return 2 * (order + 1); };
  pipeline_mng->setBoundaryRhsIntegrationRule(rule_vol);
 
  auto flux_ptr = boost::make_shared<MatrixDouble>();
  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
      new OpSetContravariantPiolaTransformOnEdge2D());
  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
      new OpCalculateHVecVectorField<3>("S", flux_ptr));
  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
      new BoundaryOp(flux_ptr, calc_flux));
 
  calc_flux = 0;
  CHKERR pipeline_mng->loopFiniteElements();
  double global_flux_assembeld = 0;
  MPI_Allreduce(&calc_flux, &global_flux_assembeld, 1, MPI_DOUBLE, MPI_SUM,
                mField.get_comm());
  calc_flux = global_flux_assembeld;
 
  MoFEMFunctionReturn(0);
}

calculateGradient()

MoFEMErrorCode Example::calculateGradient ( PetscReal *  objective_function_value, Vec  objective_function_gradient, Vec  adjoint_vector  )

private

Calculate objective function gradient using adjoint method.

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 1947 of file adjoint.cpp.

                                                              {
  MOFEM_LOG_CHANNEL("WORLD");
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  auto ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
 
  auto get_essential_fe = [this]() {
    auto post_proc_rhs = boost::make_shared<FEMethod>();
    auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
                                                           post_proc_rhs, 0)();
      MoFEMFunctionReturnHot(0);
    };
    post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
    return post_proc_rhs;
  };
 
  auto get_fd_direvative_fe = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
    // Add RHS operators for internal forces
    constexpr bool debug = false;
    if constexpr (debug) {
      auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
          mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
      pip.push_back(
          new OpAdJointTestOp<SPACE_DIM, I, DomainBaseOp>("U", common_ptr));
    } else {
      HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
          mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    }
    return fe;
  };
 
  auto calulate_fd_residual = [&](auto eps, auto diff_vec, auto fd_vec) {
    MoFEMFunctionBegin;
 
    constexpr bool debug = false;
 
    auto geom_norm = [](MoFEM::Interface &mField) {
      MoFEMFunctionBegin;
      auto field_blas = mField.getInterface<FieldBlas>();
      double nrm2 = 0.0;
      auto norm2_field = [&](const double val) {
        nrm2 += val * val;
        return val;
      };
      CHKERR field_blas->fieldLambdaOnValues(norm2_field, "GEOMETRY");
      MPI_Allreduce(MPI_IN_PLACE, &nrm2, 1, MPI_DOUBLE, MPI_SUM,
                    mField.get_comm());
      MOFEM_LOG("WORLD", Sev::inform) << "Geometry norm: " << sqrt(nrm2);
      MoFEMFunctionReturn(0);
    };
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    auto initial_current_geometry = createDMVector(adjointDM);
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
        initial_current_geometry, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecAssemblyBegin(initial_current_geometry);
    CHKERR VecAssemblyEnd(initial_current_geometry);
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    auto perturb_geometry = [&](auto eps, auto diff_vec) {
      MoFEMFunctionBegin;
      auto current_geometry = vectorDuplicate(initial_current_geometry);
      CHKERR VecCopy(initial_current_geometry, current_geometry);
      CHKERR VecAXPY(current_geometry, eps, diff_vec);
      CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
          "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
          current_geometry, INSERT_VALUES, SCATTER_REVERSE);
      MoFEMFunctionReturn(0);
    };
 
    auto fe = get_fd_direvative_fe();
    auto fp = vectorDuplicate(diff_vec);
    auto fm = vectorDuplicate(diff_vec);
    auto calc_impl = [&](auto f, auto eps) {
      MoFEMFunctionBegin;
      CHKERR VecZeroEntries(f);
      fe->f = f;
      CHKERR perturb_geometry(eps, diff_vec);
      CHKERR DMoFEMLoopFiniteElements(simple->getDM(),
                                      simple->getDomainFEName(), fe);
      CHKERR VecAssemblyBegin(f);
      CHKERR VecAssemblyEnd(f);
      CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
      auto post_proc_rhs = get_essential_fe();
      post_proc_rhs->f = f;
      CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                             post_proc_rhs.get());
      MoFEMFunctionReturn(0);
    };
    CHKERR calc_impl(fp, eps);
    CHKERR calc_impl(fm, -eps);
    CHKERR VecWAXPY(fd_vec, -1.0, fm, fp);
    CHKERR VecScale(fd_vec, 1.0 / (2.0 * eps));
 
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        "ADJOINT", "ADJOINT_FIELD", "GEOMETRY", RowColData::ROW,
        initial_current_geometry, INSERT_VALUES, SCATTER_REVERSE);
 
    if constexpr (debug)
      CHKERR geom_norm(mField);
 
    MoFEMFunctionReturn(0);
  };
 
  auto get_direvative_fe = [&](auto diff_vec) {
    auto fe_adjoint = boost::make_shared<DomainEle>(mField);
    fe_adjoint->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_adjoint->getOpPtrVector();
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
    auto cof_ptr = boost::make_shared<VectorDouble>();
 
    using OpCoFactor =
        AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
 
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "GEOMETRY", jac_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", diff_jac_ptr,
        diff_vec)); // Note: that vector is stored on displacemnt vector, that
                    // why is used here
 
    pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
 
    // Add RHS operators for internal forces
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    pip.push_back(new OpAdJointGradTimesSymTensor<SPACE_DIM, I, DomainBaseOp>(
        "U", common_ptr, jac_ptr, diff_jac_ptr, cof_ptr));
 
    return fe_adjoint;
  };
 
  auto get_objective_fe = [&](auto diff_vec, auto grad_vec,
                              auto glob_objective_ptr,
                              auto glob_objective_grad_ptr) {
    auto fe_adjoint = boost::make_shared<DomainEle>(mField);
    fe_adjoint->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data +  p_data - 1;
    };
    auto &pip = fe_adjoint->getOpPtrVector();
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    auto jac_ptr = boost::make_shared<MatrixDouble>();
    auto det_ptr = boost::make_shared<VectorDouble>();
    auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
    auto diff_jac_ptr = boost::make_shared<MatrixDouble>();
    auto cof_ptr = boost::make_shared<VectorDouble>();
    auto d_grad_ptr = boost::make_shared<MatrixDouble>();
    auto u_ptr = boost::make_shared<MatrixDouble>();
 
    using OpCoFactor =
        AdJoint<DomainEleOp>::Integration<GAUSS>::OpGetCoFactor<SPACE_DIM>;
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "GEOMETRY", jac_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", diff_jac_ptr,
        diff_vec)); // Note: that vector is stored on displacemnt vector, that
                    // why is used here
    pip.push_back(new OpCoFactor(jac_ptr, diff_jac_ptr, cof_ptr));
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "U", d_grad_ptr, grad_vec));
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, I, DomainEleOp>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
    pip.push_back(new OpAdJointObjective(
        pythonPtr, common_ptr, jac_ptr, diff_jac_ptr, cof_ptr, d_grad_ptr,
        u_ptr, glob_objective_ptr, glob_objective_grad_ptr));
 
    return fe_adjoint;
  };
 
  auto dm = simple->getDM();
  auto f = createDMVector(dm);
  auto d = vectorDuplicate(f);
  auto dm_diff_vec = vectorDuplicate(d);
 
  auto adjoint_fe = get_direvative_fe(dm_diff_vec);
  auto glob_objective_ptr = boost::make_shared<double>(0.0);
  auto glob_objective_grad_ptr = boost::make_shared<double>(0.0);
  auto objective_fe = get_objective_fe(dm_diff_vec, d, glob_objective_ptr,
                                       glob_objective_grad_ptr);
 
  auto set_variance_of_geometry = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR DMoFEMMeshToLocalVector(adjointDM, mod_vec, INSERT_VALUES,
                                   SCATTER_REVERSE);
    CHKERR mField.getInterface<VecManager>()->setOtherLocalGhostVector(
        simple->getProblemName(), "U", "ADJOINT_FIELD", RowColData::ROW,
        dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateBegin(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(dm_diff_vec, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_variance_internal_forces = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR VecZeroEntries(f);
    CHKERR VecGhostUpdateBegin(f, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(f, INSERT_VALUES, SCATTER_FORWARD);
    adjoint_fe->f = f;
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), adjoint_fe);
    CHKERR VecAssemblyBegin(f);
    CHKERR VecAssemblyEnd(f);
    CHKERR VecGhostUpdateBegin(f, ADD_VALUES, SCATTER_REVERSE);
    CHKERR VecGhostUpdateEnd(f, ADD_VALUES, SCATTER_REVERSE);
    auto post_proc_rhs = get_essential_fe();
    post_proc_rhs->f = f;
    CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                           post_proc_rhs.get());
    CHKERR VecScale(f, -1.0);
 
#ifndef NDEBUG
    constexpr bool debug = false;
    if constexpr (debug) {
      double norm0;
      CHKERR VecNorm(f, NORM_2, &norm0);
      auto fd_check = vectorDuplicate(f);
      double eps = 1e-5;
      CHKERR calulate_fd_residual(eps, dm_diff_vec, fd_check);
      double nrm;
      CHKERR VecAXPY(fd_check, -1.0, f);
      CHKERR VecNorm(fd_check, NORM_2, &nrm);
      MOFEM_LOG("WORLD", Sev::inform)
          << " FD check for internal forces [ " << mode << " ]: " << nrm
          << " / " << norm0 << " ( " << (nrm / norm0) << " )";
    }
#endif
 
    MoFEMFunctionReturn(0);
  };
 
  auto calulate_variance_of_displacement = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    CHKERR KSPSolve(kspElastic, f, d);
    CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
    CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
    MoFEMFunctionReturn(0);
  };
 
  auto calculate_variance_of_objective_function = [&](auto mode, auto mod_vec) {
    MoFEMFunctionBegin;
    *glob_objective_ptr = 0.0;
    *glob_objective_grad_ptr = 0.0;
    CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                    objective_fe);
    CHKERR VecSetValue(objective_function_gradient, mode,
                       *glob_objective_grad_ptr, ADD_VALUES);
    std::array<double, 2> array = {*glob_objective_ptr,
                                   *glob_objective_grad_ptr};
    MPI_Allreduce(MPI_IN_PLACE, array.data(), 2, MPI_DOUBLE, MPI_SUM,
                  mField.get_comm());
    *glob_objective_ptr = array[0];
    *glob_objective_grad_ptr = array[1];
    (*objective_function_value) += *glob_objective_ptr;
    MOFEM_LOG_CHANNEL("WORLD");
    MOFEM_LOG("WORLD", Sev::verbose) << " Objective gradient [ " << mode
                                     << " ]: " << *glob_objective_grad_ptr;
    MoFEMFunctionReturn(0);
  };
 
  CHKERR VecZeroEntries(objective_function_gradient);
  CHKERR VecZeroEntries(adjoint_vector);
  *objective_function_value = 0.0;
  int mode = 0;
  for (auto mod_vec : modeVecs) {
 
    CHKERR set_variance_of_geometry(mode, mod_vec);
    CHKERR calculate_variance_internal_forces(mode, mod_vec);
    CHKERR calulate_variance_of_displacement(mode, mod_vec);
    CHKERR calculate_variance_of_objective_function(mode, mod_vec);
 
    CHKERR VecAXPY(adjoint_vector, *glob_objective_grad_ptr, dm_diff_vec);
    ++mode;
  }
 
  (*objective_function_value) /= modeVecs.size();
  MOFEM_LOG("WORLD", Sev::verbose)
      << "Objective function: " << *glob_objective_ptr;
 
  CHKERR VecAssemblyBegin(objective_function_gradient);
  CHKERR VecAssemblyEnd(objective_function_gradient);
 
  CHKERR VecAssemblyBegin(adjoint_vector);
  CHKERR VecAssemblyEnd(adjoint_vector);
  CHKERR VecGhostUpdateBegin(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(adjoint_vector, INSERT_VALUES, SCATTER_FORWARD);
 
  MoFEMFunctionReturn(0);
}

checkResults() [1/10]

MoFEMErrorCode Example::checkResults ( )

private

[Postprocess results]

[Print results]

[Check]

[Check results]

[Test example]

Examples

mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-1/integration.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-8/radiation.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, and mofem/tutorials/vec-4/shallow_wave.cpp.

Definition at line 1182 of file dynamic_first_order_con_law.cpp.

                                     {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

checkResults() [2/10]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [3/10]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [4/10]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [5/10]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [6/10]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [7/10]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [8/10]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [9/10]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [10/10]

MoFEMErrorCode Example::checkResults ( )

private

createCommonData() [1/8]

MoFEMErrorCode Example::createCommonData ( )

private

[Set up problem]

[Set integration rule]

[Create common data]

< true if tau order is set

< true if tau order is set

Examples

mofem/tutorials/adv-0/plastic.cpp, mofem/tutorials/fun-1/integration.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-8/radiation.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, and mofem/tutorials/vec-4/shallow_wave.cpp.

Definition at line 480 of file plastic.cpp.

                                         {
  MoFEMFunctionBegin;
 
  auto get_command_line_parameters = [&]() {
    MoFEMFunctionBegin;
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-scale", &scale, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-young_modulus",
                                 &young_modulus, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-poisson_ratio",
                                 &poisson_ratio, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-hardening", &H, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-hardening_viscous", &visH,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-yield_stress", &sigmaY,
                                 PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn0", &cn0, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn1", &cn1, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-zeta", &zeta, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-Qinf", &Qinf, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-b_iso", &b_iso, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-C1_k", &C1_k, PETSC_NULLPTR);
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-large_strains",
                               &is_large_strains, PETSC_NULLPTR);
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-set_timer", &set_timer,
                               PETSC_NULLPTR);
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
                              PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &order, PETSC_NULLPTR);
    PetscBool tau_order_is_set; ///< true if tau order is set
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-tau_order", &tau_order,
                              &tau_order_is_set);
    PetscBool ep_order_is_set; ///< true if tau order is set
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-ep_order", &ep_order,
                              &ep_order_is_set);
    CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-geom_order", &geom_order,
                              PETSC_NULLPTR);
 
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-rho", &rho, PETSC_NULLPTR);
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-alpha_damping",
                                 &alpha_damping, PETSC_NULLPTR);
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Young modulus " << young_modulus;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Poisson ratio " << poisson_ratio;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Yield stress " << sigmaY;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Hardening " << H;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Viscous hardening " << visH;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation yield stress " << Qinf;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Saturation exponent " << b_iso;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Kinematic hardening " << C1_k;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "cn0 " << cn0;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "cn1 " << cn1;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "zeta " << zeta;
 
    if (tau_order_is_set == PETSC_FALSE)
      tau_order = order - 2;
    if (ep_order_is_set == PETSC_FALSE)
      ep_order = order - 1;
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Approximation order " << order;
    MOFEM_LOG("PLASTICITY", Sev::inform)
        << "Ep approximation order " << ep_order;
    MOFEM_LOG("PLASTICITY", Sev::inform)
        << "Tau approximation order " << tau_order;
    MOFEM_LOG("PLASTICITY", Sev::inform)
        << "Geometry approximation order " << geom_order;
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Density " << rho;
    MOFEM_LOG("PLASTICITY", Sev::inform) << "alpha_damping " << alpha_damping;
 
    PetscBool is_scale = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_scale", &is_scale,
                               PETSC_NULLPTR);
    if (is_scale) {
      scale /= young_modulus;
    }
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
 
#ifdef ADD_CONTACT
    CHKERR PetscOptionsGetScalar(PETSC_NULLPTR, "", "-cn_contact",
                                 &ContactOps::cn_contact, PETSC_NULLPTR);
    MOFEM_LOG("CONTACT", Sev::inform)
        << "cn_contact " << ContactOps::cn_contact;
#endif // ADD_CONTACT
 
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-quasi_static",
                               &is_quasi_static, PETSC_NULLPTR);
    MOFEM_LOG("PLASTICITY", Sev::inform)
        << "Is quasi static: " << (is_quasi_static ? "true" : "false");
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR get_command_line_parameters();
 
#ifdef ADD_CONTACT
  #ifdef ENABLE_PYTHON_BINDING
  auto file_exists = [](std::string myfile) {
    std::ifstream file(myfile.c_str());
    if (file) {
      return true;
    }
    return false;
  };
  char sdf_file_name[255] = "sdf.py";
  CHKERR PetscOptionsGetString(PETSC_NULLPTR, PETSC_NULLPTR, "-sdf_file",
                               sdf_file_name, 255, PETSC_NULLPTR);
 
  if (file_exists(sdf_file_name)) {
    MOFEM_LOG("CONTACT", Sev::inform) << sdf_file_name << " file found";
    sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
    CHKERR sdfPythonPtr->sdfInit(sdf_file_name);
    ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
  } else {
    MOFEM_LOG("CONTACT", Sev::warning) << sdf_file_name << " file NOT found";
  }
  #endif
#endif // ADD_CONTACT
 
  MoFEMFunctionReturn(0);
}

createCommonData() [2/8]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [3/8]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [4/8]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [5/8]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [6/8]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [7/8]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [8/8]

MoFEMErrorCode Example::createCommonData ( )

private

getCoordsInImage()

std::pair< int, int > Example::getCoordsInImage ( double  x, double  y  )

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 131 of file phase.cpp.

                                                              {
 
  auto &m = iI[focalIndex];
  x -= aveMaxMin[MIN_X];
  y -= aveMaxMin[MIN_Y];
  x *= (m.size1() - 1) / (aveMaxMin[MAX_X] - aveMaxMin[MIN_X]);
  y *= (m.size2() - 1) / (aveMaxMin[MAX_Y] - aveMaxMin[MIN_Y]);
  const auto p = std::make_pair<int, int>(std::round(x), std::round(y));
 
#ifndef NDEBUG
  if (p.first < 0 && p.first >= m.size1())
    THROW_MESSAGE("Wrong index");
  if (p.second < 0 && p.second >= m.size2())
    THROW_MESSAGE("Wrong index");
#endif
 
  return p;
}

integrateElements()

MoFEMErrorCode Example::integrateElements ( )

private

[Push operators to pipeline]

[Integrate]

Examples

mofem/tutorials/fun-1/integration.cpp.

Definition at line 218 of file integration.cpp.

                                          {
  MoFEMFunctionBegin;
  // Zero global vector
  CHKERR VecZeroEntries(commonDataPtr->petscVec);
 
  // Integrate elements by executing operators in the pipeline
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  CHKERR pipeline_mng->loopFiniteElements();
 
  // Assemble MPI vector
  CHKERR VecAssemblyBegin(commonDataPtr->petscVec);
  CHKERR VecAssemblyEnd(commonDataPtr->petscVec);
  MoFEMFunctionReturn(0);
}

integrationRule()

static int Example::integrationRule ( int  , int  , int  p_data  )

inlinestaticprivate

Examples

mofem/tutorials/scl-8/radiation.cpp.

Definition at line 52 of file radiation.cpp.

{ return 2 * p_data; };

kspSolve()

MoFEMErrorCode Example::kspSolve ( )

private

[Push operators to pipeline]

[Solve]

Examples

mofem/tutorials/scl-8/radiation.cpp.

Definition at line 230 of file radiation.cpp.

                                 {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  auto ts = pipeline_mng->createTSIM();
 
  double ftime = 1;
  CHKERR TSSetDuration(ts, PETSC_DEFAULT, ftime);
  CHKERR TSSetFromOptions(ts);
  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
 
  auto T = createDMVector(simple->getDM());
  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                 SCATTER_FORWARD);
 
  CHKERR TSSolve(ts, T);
  CHKERR TSGetTime(ts, &ftime);
 
  PetscInt steps, snesfails, rejects, nonlinits, linits;
  CHKERR TSGetTimeStepNumber(ts, &steps);
  CHKERR TSGetSNESFailures(ts, &snesfails);
  CHKERR TSGetStepRejections(ts, &rejects);
  CHKERR TSGetSNESIterations(ts, &nonlinits);
  CHKERR TSGetKSPIterations(ts, &linits);
  MOFEM_LOG_C("EXAMPLE", Sev::inform,
              "steps %d (%d rejected, %d SNES fails), ftime %g, nonlinits "
              "%d, linits %d",
              steps, rejects, snesfails, ftime, nonlinits, linits);
 
  MoFEMFunctionReturn(0);
}

lhsFlux()

double Example::lhsFlux ( const double  x, const double  y, const double    )

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 165 of file phase.cpp.

                                                                    {
  const auto idx = getCoordsInImage(x, y);
  const auto &m = iI[focalIndex];
  return 1. / m(idx.first, idx.second);
}

OPs() [1/2]

MoFEMErrorCode Example::OPs ( )

private

[Boundary condition]

[Push operators to pipeline]

[Only used for dynamics]

[Only used for dynamics]

[Only used for dynamics]

[Only used for dynamics]

Examples

mofem/tutorials/adv-0/plastic.cpp, and mofem/tutorials/scl-8/radiation.cpp.

Definition at line 650 of file plastic.cpp.

                            {
  MoFEMFunctionBegin;
  auto pip_mng = mField.getInterface<PipelineManager>();
 
  auto integration_rule_bc = [](int, int, int ao) { return 2 * ao; };
 
  auto vol_rule = [](int, int, int ao) { return 2 * ao + geom_order - 1; };
 
  auto add_boundary_ops_lhs_mechanical = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR PlasticOps::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        pip, {HDIV}, "GEOMETRY");
    pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
 
    // Add Natural BCs to LHS
    CHKERR BoundaryLhsBCs::AddFluxToPipeline<OpBoundaryLhsBCs>::add(
        pip, mField, "U", Sev::inform);
 
#ifdef ADD_CONTACT
    auto simple = mField.getInterface<Simple>();
    CHKERR
    ContactOps::opFactoryBoundaryLhs<SPACE_DIM, AT, GAUSS, BoundaryEleOp>(
        pip, "SIGMA", "U");
    CHKERR
    ContactOps::opFactoryBoundaryToDomainLhs<SPACE_DIM, AT, IT, DomainEle>(
        mField, pip, simple->getDomainFEName(), "SIGMA", "U", "GEOMETRY",
        vol_rule);
#endif // ADD_CONTACT
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_boundary_ops_rhs_mechanical = [&](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR PlasticOps::AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        pip, {HDIV}, "GEOMETRY");
    pip.push_back(new OpSetHOWeightsOnSubDim<SPACE_DIM>());
 
    // Add Natural BCs to RHS
    CHKERR BoundaryRhsBCs::AddFluxToPipeline<OpBoundaryRhsBCs>::add(
        pip, mField, "U", {boost::make_shared<ScaledTimeScale>()}, Sev::inform);
 
#ifdef ADD_CONTACT
    CHKERR ContactOps::opFactoryBoundaryRhs<SPACE_DIM, AT, IT, BoundaryEleOp>(
        pip, "SIGMA", "U");
#endif // ADD_CONTACT
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_lhs = [this](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        pip, {H1, HDIV}, "GEOMETRY");
 
    if (is_quasi_static == PETSC_FALSE) {
 
      //! [Only used for dynamics]
      using OpMass = FormsIntegrators<DomainEleOp>::Assembly<AT>::BiLinearForm<
          GAUSS>::OpMass<1, SPACE_DIM>;
      //! [Only used for dynamics]
 
      auto get_inertia_and_mass_damping = [this](const double, const double,
                                                 const double) {
        auto *pip = mField.getInterface<PipelineManager>();
        auto &fe_domain_lhs = pip->getDomainLhsFE();
        return (rho / scale) * fe_domain_lhs->ts_aa +
               (alpha_damping / scale) * fe_domain_lhs->ts_a;
      };
      pip.push_back(new OpMass("U", "U", get_inertia_and_mass_damping));
    }
 
    CHKERR PlasticOps::opFactoryDomainLhs<SPACE_DIM, AT, IT, DomainEleOp>(
        mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_domain_ops_rhs = [this](auto &pip) {
    MoFEMFunctionBegin;
 
    CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        pip, {H1, HDIV}, "GEOMETRY");
 
    CHKERR DomainRhsBCs::AddFluxToPipeline<OpDomainRhsBCs>::add(
        pip, mField, "U",
        {boost::make_shared<ScaledTimeScale>("body_force_hist.txt")},
        Sev::inform);
 
    // only in case of dynamics
    if (is_quasi_static == PETSC_FALSE) {
 
      //! [Only used for dynamics]
      using OpInertiaForce = FormsIntegrators<DomainEleOp>::Assembly<
          AT>::LinearForm<IT>::OpBaseTimesVector<1, SPACE_DIM, 1>;
      //! [Only used for dynamics]
 
      auto mat_acceleration = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValuesDotDot<SPACE_DIM>(
          "U", mat_acceleration));
      pip.push_back(
          new OpInertiaForce("U", mat_acceleration, [](double, double, double) {
            return rho / scale;
          }));
      if (alpha_damping > 0) {
        auto mat_velocity = boost::make_shared<MatrixDouble>();
        pip.push_back(
            new OpCalculateVectorFieldValuesDot<SPACE_DIM>("U", mat_velocity));
        pip.push_back(
            new OpInertiaForce("U", mat_velocity, [](double, double, double) {
              return alpha_damping / scale;
            }));
      }
    }
 
    CHKERR PlasticOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
        mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
 
#ifdef ADD_CONTACT
    CHKERR ContactOps::opFactoryDomainRhs<SPACE_DIM, AT, IT, DomainEleOp>(
        pip, "SIGMA", "U");
#endif // ADD_CONTACT
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR add_domain_ops_lhs(pip_mng->getOpDomainLhsPipeline());
  CHKERR add_domain_ops_rhs(pip_mng->getOpDomainRhsPipeline());
 
  // Boundary
  CHKERR add_boundary_ops_lhs_mechanical(pip_mng->getOpBoundaryLhsPipeline());
  CHKERR add_boundary_ops_rhs_mechanical(pip_mng->getOpBoundaryRhsPipeline());
 
  CHKERR pip_mng->setDomainRhsIntegrationRule(vol_rule);
  CHKERR pip_mng->setDomainLhsIntegrationRule(vol_rule);
 
  CHKERR pip_mng->setBoundaryLhsIntegrationRule(integration_rule_bc);
  CHKERR pip_mng->setBoundaryRhsIntegrationRule(integration_rule_bc);
 
  auto create_reaction_pipeline = [&](auto &pip) {
    MoFEMFunctionBegin;
    CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        pip, {H1}, "GEOMETRY");
    CHKERR PlasticOps::opFactoryDomainReactions<SPACE_DIM, AT, IT, DomainEleOp>(
        mField, "MAT_PLASTIC", pip, "U", "EP", "TAU");
    MoFEMFunctionReturn(0);
  };
 
  reactionFe = boost::make_shared<DomainEle>(mField);
  reactionFe->getRuleHook = vol_rule;
  CHKERR create_reaction_pipeline(reactionFe->getOpPtrVector());
  reactionFe->postProcessHook =
      EssentialPreProcReaction<DisplacementCubitBcData>(mField, reactionFe);
 
  MoFEMFunctionReturn(0);
}

OPs() [2/2]

MoFEMErrorCode Example::OPs ( )

private

outputResults() [1/9]

MoFEMErrorCode Example::outputResults ( )

private

[Solve]

[Postprocess results]

Examples

mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/mix-1/phase.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, and mofem/tutorials/vec-4/shallow_wave.cpp.

Definition at line 1160 of file dynamic_first_order_con_law.cpp.

                                      {
  MoFEMFunctionBegin;
  PetscBool test_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test", &test_flg, PETSC_NULLPTR);
  if (test_flg) {
    auto *simple = mField.getInterface<Simple>();
    auto T = createDMVector(simple->getDM());
    CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                   SCATTER_FORWARD);
    double nrm2;
    CHKERR VecNorm(T, NORM_2, &nrm2);
    MOFEM_LOG("EXAMPLE", Sev::inform) << "Regression norm " << nrm2;
    constexpr double regression_value = 0.0194561;
    if (fabs(nrm2 - regression_value) > 1e-2)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
              "Regression test failed; wrong norm value.");
  }
  MoFEMFunctionReturn(0);
}

outputResults() [2/9]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [3/9]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [4/9]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [5/9]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [6/9]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [7/9]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [8/9]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [9/9]

MoFEMErrorCode Example::outputResults ( const int  i )

private

[Solve]

[Postprocess results]

Definition at line 504 of file phase.cpp.

                                                 {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
  pipeline_mng->getDomainLhsFE().reset();
  pipeline_mng->getDomainRhsFE().reset();
  pipeline_mng->getBoundaryRhsFE().reset();
 
  using PostProcEle = PostProcBrokenMeshInMoab<DomainEle>;
 
  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateHOJacForFace(jac_ptr));
  post_proc_fe->getOpPtrVector().push_back(new OpMakeHdivFromHcurl());
  post_proc_fe->getOpPtrVector().push_back(
      new OpSetContravariantPiolaTransformOnFace2D(jac_ptr));
 
  auto s_ptr = boost::make_shared<VectorDouble>();
  auto phi_ptr = boost::make_shared<MatrixDouble>();
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateScalarFieldValues("S", s_ptr));
  post_proc_fe->getOpPtrVector().push_back(
      new OpCalculateHVecVectorField<3>("PHI", phi_ptr));
 
  using OpPPMap = OpPostProcMapInMoab<3, 3>;
 
  post_proc_fe->getOpPtrVector().push_back(
 
      new OpPPMap(post_proc_fe->getPostProcMesh(),
                  post_proc_fe->getMapGaussPts(),
 
                  OpPPMap::DataMapVec{{"S", s_ptr}},
 
                  OpPPMap::DataMapMat{{"PHI", phi_ptr}},
 
                  OpPPMap::DataMapMat{},
 
                  OpPPMap::DataMapMat{}
 
                  )
 
  );
 
  // post_proc_fe->addFieldValuesPostProc("S");
  // post_proc_fe->addFieldValuesPostProc("PHI");
 
 
  pipeline_mng->getDomainRhsFE() = post_proc_fe;
  CHKERR pipeline_mng->loopFiniteElements();
  CHKERR post_proc_fe->writeFile("out_" + boost::lexical_cast<std::string>(i) +
                                 ".h5m");
  MoFEMFunctionReturn(0);
}

postProcess() [1/2]

MoFEMErrorCode Example::postProcess ( )

private

[Integrate]

[Solve]

[Print results]

[Postprocess results]

Examples

mofem/tutorials/fun-1/integration.cpp, and mofem/tutorials/scl-8/radiation.cpp.

Definition at line 235 of file integration.cpp.

                                    {
  MoFEMFunctionBegin;
  const double *array;
  CHKERR VecGetArrayRead(commonDataPtr->petscVec, &array);
  if (mField.get_comm_rank() == 0) {
    MOFEM_LOG_C("SELF", Sev::inform, "Mass %6.4e", array[CommonData::ZERO]);
    MOFEM_LOG_C("SELF", Sev::inform,
                "First moment of inertia [ %6.4e, %6.4e, %6.4e ]",
                array[CommonData::FIRST_X], array[CommonData::FIRST_Y],
                array[CommonData::FIRST_Z]);
    MOFEM_LOG_C("SELF", Sev::inform,
                "Second moment of inertia [ %6.4e, %6.4e, %6.4e; %6.4e %6.4e; "
                "%6.4e ]",
                array[CommonData::SECOND_XX], array[CommonData::SECOND_XY],
                array[CommonData::SECOND_XZ], array[CommonData::SECOND_YY],
                array[CommonData::SECOND_YZ], array[CommonData::SECOND_ZZ]);
  }
  CHKERR VecRestoreArrayRead(commonDataPtr->petscVec, &array);
  MoFEMFunctionReturn(0);
}

postProcess() [2/2]

MoFEMErrorCode Example::postProcess ( )

private

postprocessElastic()

MoFEMErrorCode Example::postprocessElastic ( int  iter, SmartPetscObj< Vec >  adjoint_vector = nullptr  )

private

Post-process and output results.

[Solve]

[Postprocess results]

[Postprocess clean]

[Postprocess clean]

[Postprocess initialise]

[Postprocess initialise]

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 1347 of file adjoint.cpp.

                                                                              {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto pip = mField.getInterface<PipelineManager>();
  auto det_ptr = boost::make_shared<VectorDouble>();
  auto jac_ptr = boost::make_shared<MatrixDouble>();
  auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
  //! [Postprocess clean]
  pip->getDomainRhsFE().reset();
  pip->getDomainLhsFE().reset();
  pip->getBoundaryRhsFE().reset();
  pip->getBoundaryLhsFE().reset();
  //! [Postprocess clean]
 
  //! [Postprocess initialise]
  auto post_proc_mesh = boost::make_shared<moab::Core>();
  auto post_proc_begin = boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(
      mField, post_proc_mesh);
  auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
      mField, post_proc_mesh);
  //! [Postprocess initialise]
 
  // lambada function to calculate displacement, strain and stress
  auto calculate_stress_ops = [&](auto &pip) {
    // Add H1 geometry ops
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1}, "GEOMETRY");
 
    // Use HookeOps commonDataFactory to get matStrainPtr and matCauchyStressPtr
    auto common_ptr = HookeOps::commonDataFactory<SPACE_DIM, GAUSS, DomainEle>(
        mField, pip, "U", "MAT_ELASTIC", Sev::noisy);
 
    // Store U and GEOMETRY values if needed
    auto u_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
    auto x_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
 
    using DataMapMat = std::map<std::string, boost::shared_ptr<MatrixDouble>>;
 
    DataMapMat vec_map = {{"U", u_ptr}, {"GEOMETRY", x_ptr}};
    DataMapMat sym_map = {{"STRAIN", common_ptr->getMatStrain()},
                          {"STRESS", common_ptr->getMatCauchyStress()}};
 
    if (adjoint_vector) {
      auto adjoint_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>(
          "U", adjoint_ptr, adjoint_vector));
      vec_map["ADJOINT"] = adjoint_ptr;
    }
 
    // Return what you need: displacements, coordinates, strain, stress
    return boost::make_tuple(u_ptr, x_ptr, common_ptr->getMatStrain(),
                             common_ptr->getMatCauchyStress(), vec_map,
                             sym_map);
  };
 
  auto get_tag_id_on_pmesh = [&](bool post_proc_skin) {
    int def_val_int = 0;
    Tag tag_mat;
    CHKERR mField.get_moab().tag_get_handle(
        "MAT_ELASTIC", 1, MB_TYPE_INTEGER, tag_mat,
        MB_TAG_CREAT | MB_TAG_SPARSE, &def_val_int);
    auto meshset_vec_ptr =
        mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(
            std::regex((boost::format("%s(.*)") % "MAT_ELASTIC").str()));
 
    Range skin_ents;
    std::unique_ptr<Skinner> skin_ptr;
    if (post_proc_skin) {
      skin_ptr = std::make_unique<Skinner>(&mField.get_moab());
      auto boundary_meshset = simple->getBoundaryMeshSet();
      CHKERR mField.get_moab().get_entities_by_handle(boundary_meshset,
                                                      skin_ents, true);
    }
 
    for (auto m : meshset_vec_ptr) {
      Range ents_3d;
      CHKERR mField.get_moab().get_entities_by_handle(m->getMeshset(), ents_3d,
                                                      true);
      int const id = m->getMeshsetId();
      ents_3d = ents_3d.subset_by_dimension(SPACE_DIM);
      if (post_proc_skin) {
        Range skin_faces;
        CHKERR skin_ptr->find_skin(0, ents_3d, false, skin_faces);
        ents_3d = intersect(skin_ents, skin_faces);
      }
      CHKERR mField.get_moab().tag_clear_data(tag_mat, ents_3d, &id);
    }
 
    return tag_mat;
  };
 
  auto post_proc_domain = [&](auto post_proc_mesh) {
    auto post_proc_fe =
        boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
        calculate_stress_ops(post_proc_fe->getOpPtrVector());
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            vec_map,
 
            {},
 
            sym_map
 
            )
 
    );
 
    post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(false)});
    return post_proc_fe;
  };
 
  auto post_proc_boundary = [&](auto post_proc_mesh) {
    auto post_proc_fe =
        boost::make_shared<PostProcEleBdy>(mField, post_proc_mesh);
    AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
    auto op_loop_side =
        new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
    // push ops to side element, through op_loop_side operator
    auto [u_ptr, x_ptr, mat_strain_ptr, mat_stress_ptr, vec_map, sym_map] =
        calculate_stress_ops(op_loop_side->getOpPtrVector());
    post_proc_fe->getOpPtrVector().push_back(op_loop_side);
    auto mat_traction_ptr = boost::make_shared<MatrixDouble>();
    post_proc_fe->getOpPtrVector().push_back(
        new ElasticOps::OpCalculateTraction(mat_stress_ptr, mat_traction_ptr));
    vec_map["T"] = mat_traction_ptr;
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            vec_map,
 
            {},
 
            sym_map
 
            )
 
    );
 
    post_proc_fe->setTagsToTransfer({get_tag_id_on_pmesh(true)});
    return post_proc_fe;
  };
 
  PetscBool post_proc_skin_only = PETSC_FALSE;
  if (SPACE_DIM == 3) {
    post_proc_skin_only = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin_only",
                               &post_proc_skin_only, PETSC_NULLPTR);
  }
  if (post_proc_skin_only == PETSC_FALSE) {
    pip->getDomainRhsFE() = post_proc_domain(post_proc_mesh);
  } else {
    pip->getBoundaryRhsFE() = post_proc_boundary(post_proc_mesh);
  }
 
  CHKERR DMoFEMPreProcessFiniteElements(simple->getDM(),
                                        post_proc_begin->getFEMethod());
  CHKERR pip->loopFiniteElements();
  CHKERR DMoFEMPostProcessFiniteElements(simple->getDM(),
                                         post_proc_end->getFEMethod());
 
  CHKERR post_proc_end->writeFile("out_elastic_" + std::to_string(iter) +
                                  ".h5m");
  MoFEMFunctionReturn(0);
}

pushOperators()

MoFEMErrorCode Example::pushOperators ( )

private

[Set density distribution]

[Push operators to pipeline]

Examples

mofem/tutorials/fun-1/integration.cpp.

Definition at line 188 of file integration.cpp.

                                      {
  MoFEMFunctionBegin;
  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
 
  // Push an operator which calculates values of density at integration points
  pipeline_mng->getOpDomainRhsPipeline().push_back(
      new OpCalculateScalarFieldValues(
          "rho", commonDataPtr->getRhoAtIntegrationPtsPtr()));
 
  // Push an operator to pipeline to calculate zero moment of inertia (mass)
  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpZero(commonDataPtr));
 
  // Push an operator to the pipeline to calculate first moment of inertaia
  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpFirst(commonDataPtr));
 
  // Push an operator to the pipeline to calculate second moment of inertaia
  pipeline_mng->getOpDomainRhsPipeline().push_back(new OpSecond(commonDataPtr));
 
  // Set integration rule. Integration rule is equal to the polynomial order of
  // the density field plus 2, since under the integral of the second moment of
  // inertia term x*x is present
  auto integration_rule = [](int, int, int p_data) { return p_data + 2; };
 
  // Add integration rule to the element
  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
  MoFEMFunctionReturn(0);
}

readMesh() [1/10]

MoFEMErrorCode Example::readMesh ( )

private

[Run problem]

[Run programme]

[run problem]

[Read mesh]

[Read mesh]

Read mesh from file and setup material/boundary condition meshsets

This function loads the finite element mesh from file and processes associated meshsets that define material properties and boundary conditions. The mesh is typically generated using CUBIT and exported in .h5m format.

Meshsets are used to group elements/faces by:

• Material properties (for different material blocks)
• Boundary conditions (for applying loads and constraints)
• Optimization regions (for topology optimization)

Returns

MoFEMErrorCode Success or error code

Examples

mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/mix-1/phase.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4/shallow_wave.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 462 of file dynamic_first_order_con_law.cpp.

                                 {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
  MoFEMFunctionReturn(0);
}

readMesh() [2/10]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [3/10]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [4/10]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [5/10]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [6/10]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [7/10]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [8/10]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [9/10]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [10/10]

MoFEMErrorCode Example::readMesh ( )

private

Read mesh from file and setup meshsets.

rhsSource()

double Example::rhsSource ( const double  x, const double  y, const double    )

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 150 of file phase.cpp.

                                                                      {
  const auto idx = getCoordsInImage(x, y);
 
  double v = 0;
  for (auto w = 0; w != window_savitzky_golay; ++w) {
    const auto i = focalIndex - (window_savitzky_golay - 1) / 2 + w;
    const auto &intensity = iI[i];
    v += intensity(idx.first, idx.second) * savitzkyGolayWeights[w];
  }
  v = static_cast<double>(v) / savitzkyGolayNormalisation;
 
  const auto dz = rZ[focalIndex + 1] - rZ[focalIndex - 1];
  return -k * v / dz;
}

runProblem() [1/13]

MoFEMErrorCode Example::runProblem

(

)

[Run problem]

[Run topology optimization problem]

[Create common data]

[Run programme]

[Operator]

[run problem]

[Run all]

[Run problem]

Main driver for topology optimization using adjoint sensitivity analysis

This function orchestrates the complete topology optimization workflow:

1. Initialize Python objective function interface
2. Setup finite element problems (forward and adjoint)
3. Compute initial elastic solution
4. Generate topology optimization modes
5. Run TAO optimization loop with adjoint-based gradients

The optimization uses TAO (Toolkit for Advanced Optimization) with L-BFGS algorithm. At each iteration:

• Update geometry based on current design variables
• Solve forward elastic problem
  
• Compute objective function and gradients using adjoint method
• Post-process results

Returns

MoFEMErrorCode Success or error code

Setup TAO (Toolkit for Advanced Optimization) solver for topology optimization

TAO provides various optimization algorithms. Here we use TAOLMVM (Limited Memory Variable Metric) which is a quasi-Newton method suitable for unconstrained optimization with gradient information.

The optimization variables are coefficients for the topology modes, and gradients are computed using the adjoint method for efficiency.

Generate modes for topology optimization design parameterization

These modes represent perturbations to the geometry that can be used as design variables in topology optimization. The modes are defined through the Python interface and provide spatial basis functions for shape modifications.

Start optimization with zero initial guess for design variables

The TAO solver will iteratively:

1. Evaluate objective function at current design point
2. Compute gradients using adjoint sensitivity analysis
   
3. Update design variables using L-BFGS algorithm
4. Check convergence criteria

Each function evaluation involves solving the forward elastic problem and the adjoint problem to compute sensitivities efficiently.

Examples

mofem/tutorials/adv-0/plastic.cpp, mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-1/integration.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/mix-1/phase.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-8/radiation.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4/shallow_wave.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 256 of file plastic.cpp.

                                   {
  MoFEMFunctionBegin;
  CHKERR createCommonData();
  CHKERR setupProblem();
  CHKERR bC();
  CHKERR OPs();
  PetscBool test_ops = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-test_operators", &test_ops,
                             PETSC_NULLPTR);
  if (test_ops == PETSC_FALSE) {
    CHKERR tsSolve();
  } else {
    CHKERR testOperators();
  }
  MoFEMFunctionReturn(0);
}

runProblem() [2/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [3/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [4/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [5/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [6/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [7/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [8/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [9/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [10/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [11/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [12/13]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [13/13]

MoFEMErrorCode Example::runProblem

(

)

Main driver function for the optimization process.

setFieldValues()

MoFEMErrorCode Example::setFieldValues ( )

private

[Create common data]

[Set density distribution]

Examples

mofem/tutorials/fun-1/integration.cpp.

Definition at line 172 of file integration.cpp.

                                       {
  MoFEMFunctionBegin;
  auto set_density = [&](VectorAdaptor &&field_data, double *xcoord,
                         double *ycoord, double *zcoord) {
    MoFEMFunctionBeginHot;
    field_data[0] = 1;
    MoFEMFunctionReturnHot(0);
  };
  FieldBlas *field_blas;
  CHKERR mField.getInterface(field_blas);
  CHKERR field_blas->setVertexDofs(set_density, "rho");
  MoFEMFunctionReturn(0);
}

setIntegrationRules() [1/3]

MoFEMErrorCode Example::setIntegrationRules ( )

private

[Set up problem]

[Set integration rule]

Examples

mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/scl-0/approximaton.cpp, and mofem/tutorials/scl-9/heat_method.cpp.

Definition at line 229 of file plot_base.cpp.

                                            {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

setIntegrationRules() [2/3]

MoFEMErrorCode Example::setIntegrationRules ( )

private

setIntegrationRules() [3/3]

MoFEMErrorCode Example::setIntegrationRules ( )

private

setUp()

MoFEMErrorCode Example::setUp ( )

private

[Run all]

[Set up problem]

Examples

mofem/tutorials/fun-1/integration.cpp.

Definition at line 137 of file integration.cpp.

                              {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
  CHKERR simple->getOptions();
  CHKERR simple->loadFile();
  // Add field
  CHKERR simple->addDomainField("rho", H1, AINSWORTH_LEGENDRE_BASE, 1);
  constexpr int order = 1;
  CHKERR simple->setFieldOrder("rho", order);
  CHKERR simple->setUp();
  MoFEMFunctionReturn(0);
}

setupAdJoint()

MoFEMErrorCode Example::setupAdJoint ( )

private

Setup adjoint fields and finite elements

[Set up problem]

[Setup adjoint]

Setup adjoint fields and finite elements for sensitivity analysis

The adjoint method is used to efficiently compute gradients of the objective function with respect to design variables. This function sets up:

1. ADJOINT_FIELD - stores adjoint variables (Lagrange multipliers)
2. ADJOINT_DM - data manager for adjoint problem
3. Adjoint finite elements for domain and boundary

The adjoint equation is: K^T * λ = ∂f/∂u where λ are adjoint variables, K is stiffness matrix, f is objective

Returns

MoFEMErrorCode Success or error code

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 623 of file adjoint.cpp.

                                     {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
 
  // Create adjoint data manager and field
  auto create_adjoint_dm = [&]() {
    auto adjoint_dm = createDM(mField.get_comm(), "DMMOFEM");
 
    auto add_field = [&]() {
      MoFEMFunctionBegin;
      CHKERR mField.add_field("ADJOINT_FIELD", H1, base, SPACE_DIM);
      CHKERR mField.add_ents_to_field_by_dim(simple->getMeshset(), SPACE_DIM,
                                             "ADJOINT_FIELD");
      for (auto tt = MBEDGE; tt <= moab::CN::TypeDimensionMap[SPACE_DIM].second;
           ++tt)
        CHKERR mField.set_field_order(simple->getMeshset(), tt, "ADJOINT_FIELD",
                                      fieldOrder);
      CHKERR mField.set_field_order(simple->getMeshset(), MBVERTEX,
                                    "ADJOINT_FIELD", 1);
      CHKERR mField.build_fields();
      MoFEMFunctionReturn(0);
    };
 
    auto add_adjoint_fe_impl = [&]() {
      MoFEMFunctionBegin;
      CHKERR mField.add_finite_element("ADJOINT_DOMAIN_FE");
      CHKERR mField.modify_finite_element_add_field_row("ADJOINT_DOMAIN_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_col("ADJOINT_DOMAIN_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_DOMAIN_FE",
                                                         "GEOMETRY");
      CHKERR mField.add_ents_to_finite_element_by_dim(
          simple->getMeshset(), SPACE_DIM, "ADJOINT_DOMAIN_FE");
      CHKERR mField.build_finite_elements("ADJOINT_DOMAIN_FE");
 
      CHKERR mField.add_finite_element("ADJOINT_BOUNDARY_FE");
      CHKERR mField.modify_finite_element_add_field_row("ADJOINT_BOUNDARY_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_col("ADJOINT_BOUNDARY_FE",
                                                        "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "ADJOINT_FIELD");
      CHKERR mField.modify_finite_element_add_field_data("ADJOINT_BOUNDARY_FE",
                                                         "GEOMETRY");
 
      auto block_name = "OPTIMISE";
      auto mesh_mng = mField.getInterface<MeshsetsManager>();
      auto bcs = mesh_mng->getCubitMeshsetPtr(
 
          std::regex((boost::format("%s(.*)") % block_name).str())
 
      );
 
      for (auto bc : bcs) {
        CHKERR mField.add_ents_to_finite_element_by_dim(
            bc->getMeshset(), SPACE_DIM - 1, "ADJOINT_BOUNDARY_FE");
      }
 
      CHKERR mField.build_finite_elements("ADJOINT_BOUNDARY_FE");
 
      CHKERR mField.build_adjacencies(simple->getBitRefLevel(),
                                      simple->getBitRefLevelMask());
 
      MoFEMFunctionReturn(0);
    };
 
    auto set_adjoint_dm_imp = [&]() {
      MoFEMFunctionBegin;
      CHKERR DMMoFEMCreateMoFEM(adjoint_dm, &mField, "ADJOINT",
                                simple->getBitRefLevel(),
                                simple->getBitRefLevelMask());
      CHKERR DMMoFEMSetDestroyProblem(adjoint_dm, PETSC_TRUE);
      CHKERR DMSetFromOptions(adjoint_dm);
      CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_DOMAIN_FE");
      CHKERR DMMoFEMAddElement(adjoint_dm, "ADJOINT_BOUNDARY_FE");
      CHKERR DMMoFEMSetSquareProblem(adjoint_dm, PETSC_TRUE);
      CHKERR DMMoFEMSetIsPartitioned(adjoint_dm, PETSC_TRUE);
      mField.getInterface<ProblemsManager>()->buildProblemFromFields =
          PETSC_TRUE;
      CHKERR DMSetUp(adjoint_dm);
      mField.getInterface<ProblemsManager>()->buildProblemFromFields =
          PETSC_FALSE;
      MoFEMFunctionReturn(0);
    };
 
    CHK_THROW_MESSAGE(add_field(), "add adjoint field");
    CHK_THROW_MESSAGE(add_adjoint_fe_impl(), "add adjoint fe");
    CHK_THROW_MESSAGE(set_adjoint_dm_imp(), "set adjoint dm");
 
    return adjoint_dm;
  };
 
  adjointDM = create_adjoint_dm();
 
  MoFEMFunctionReturn(0);
}

setupProblem() [1/12]

MoFEMErrorCode Example::setupProblem ( )

private

[Run problem]

[Read mesh]

[Set up problem]

[Set up problem]

Setup finite element fields, approximation spaces and degrees of freedom

This function configures the finite element problem by:

1. Setting up the displacement field "U" with vector approximation
2. Setting up the geometry field "GEOMETRY" for mesh deformation
   
3. Defining polynomial approximation order and basis functions
4. Creating degrees of freedom on mesh entities

The displacement field uses H1 vector space for standard elasticity. The geometry field allows mesh modification during topology optimization. Different basis functions (Ainsworth-Legendre vs Demkowicz) can be selected.

Returns

MoFEMErrorCode Success or error code

Setup displacement field "U" - the primary unknown in elasticity This field represents displacement vector at each node/DOF

Setup geometry field "GEOMETRY" - used for mesh deformation in optimization This field stores current nodal coordinates and can be modified during topology optimization to represent design changes

For higher-order elements, this projects the exact geometry onto the geometry field to maintain curved boundaries accurately

Examples

mofem/tutorials/adv-0/plastic.cpp, mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/mix-1/phase.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-8/radiation.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4/shallow_wave.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 275 of file plastic.cpp.

                                     {
  MoFEMFunctionBegin;
  Simple *simple = mField.getInterface<Simple>();
 
  Range domain_ents;
  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
                                                     true);
  auto get_ents_by_dim = [&](const auto dim) {
    if (dim == SPACE_DIM) {
      return domain_ents;
    } else {
      Range ents;
      if (dim == 0)
        CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
      else
        CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
      return ents;
    }
  };
 
  auto get_base = [&]() {
    auto domain_ents = get_ents_by_dim(SPACE_DIM);
    if (domain_ents.empty())
      CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
    const auto type = type_from_handle(domain_ents[0]);
    switch (type) {
    case MBQUAD:
      return DEMKOWICZ_JACOBI_BASE;
    case MBHEX:
      return DEMKOWICZ_JACOBI_BASE;
    case MBTRI:
      return AINSWORTH_LEGENDRE_BASE;
    case MBTET:
      return AINSWORTH_LEGENDRE_BASE;
    default:
      CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
    }
    return NOBASE;
  };
 
  const auto base = get_base();
  MOFEM_LOG("PLASTICITY", Sev::inform)
      << "Base " << ApproximationBaseNames[base];
 
  CHKERR simple->addDomainField("U", H1, base, SPACE_DIM);
  CHKERR simple->addDomainField("EP", L2, base, size_symm);
  CHKERR simple->addDomainField("TAU", L2, base, 1);
  CHKERR simple->addBoundaryField("U", H1, base, SPACE_DIM);
 
  CHKERR simple->addDataField("GEOMETRY", H1, base, SPACE_DIM);
 
  PetscBool order_edge = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_edge", &order_edge,
                             PETSC_NULLPTR);
  PetscBool order_face = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_face", &order_face,
                             PETSC_NULLPTR);
  PetscBool order_volume = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-order_volume", &order_volume,
                             PETSC_NULLPTR);
 
  if (order_edge || order_face || order_volume) {
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Order edge " << order_edge
        ? "true"
        : "false";
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Order face " << order_face
        ? "true"
        : "false";
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Order volume " << order_volume
        ? "true"
        : "false";
 
    auto ents = get_ents_by_dim(0);
    if (order_edge)
      ents.merge(get_ents_by_dim(1));
    if (order_face)
      ents.merge(get_ents_by_dim(2));
    if (order_volume)
      ents.merge(get_ents_by_dim(3));
    CHKERR simple->setFieldOrder("U", order, &ents);
  } else {
    CHKERR simple->setFieldOrder("U", order);
  }
  CHKERR simple->setFieldOrder("EP", ep_order);
  CHKERR simple->setFieldOrder("TAU", tau_order);
 
  CHKERR simple->setFieldOrder("GEOMETRY", geom_order);
 
#ifdef ADD_CONTACT
  CHKERR simple->addDomainField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                SPACE_DIM);
  CHKERR simple->addBoundaryField("SIGMA", CONTACT_SPACE, DEMKOWICZ_JACOBI_BASE,
                                  SPACE_DIM);
 
  auto get_skin = [&]() {
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents);
    Skinner skin(&mField.get_moab());
    Range skin_ents;
    CHKERR skin.find_skin(0, body_ents, false, skin_ents);
    return skin_ents;
  };
 
  auto filter_blocks = [&](auto skin) {
    bool is_contact_block = true;
    Range contact_range;
    for (auto m :
         mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
             (boost::format("%s(.*)") % "CONTACT").str()
 
                 ))
 
    ) {
      is_contact_block =
          true; ///< blocs interation is collective, so that is set irrespective
                ///< if there are entities in given rank or not in the block
      MOFEM_LOG("CONTACT", Sev::inform)
          << "Find contact block set:  " << m->getName();
      auto meshset = m->getMeshset();
      Range contact_meshset_range;
      CHKERR mField.get_moab().get_entities_by_dimension(
          meshset, SPACE_DIM - 1, contact_meshset_range, true);
 
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          contact_meshset_range);
      contact_range.merge(contact_meshset_range);
    }
    if (is_contact_block) {
      MOFEM_LOG("SYNC", Sev::inform)
          << "Nb entities in contact surface: " << contact_range.size();
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
      skin = intersect(skin, contact_range);
    }
    return skin;
  };
 
  auto filter_true_skin = [&](auto skin) {
    Range boundary_ents;
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(skin, PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, &boundary_ents);
    return boundary_ents;
  };
 
  auto boundary_ents = filter_true_skin(filter_blocks(get_skin()));
  CHKERR simple->setFieldOrder("SIGMA", 0);
  CHKERR simple->setFieldOrder("SIGMA", order - 1, &boundary_ents);
#endif
 
  CHKERR simple->setUp();
  CHKERR simple->addFieldToEmptyFieldBlocks("U", "TAU");
 
  auto project_ho_geometry = [&]() {
    Projection10NodeCoordsOnField ent_method(mField, "GEOMETRY");
    return mField.loop_dofs("GEOMETRY", ent_method);
  };
  PetscBool project_geometry = PETSC_TRUE;
  CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-project_geometry",
                             &project_geometry, PETSC_NULLPTR);
  if (project_geometry) {
    CHKERR project_ho_geometry();
  }
 
  auto get_volume = [&]() {
    using VolOp = DomainEle::UserDataOperator;
    auto *op_ptr = new VolOp(NOSPACE, VolOp::OPSPACE);
    std::array<double, 2> volume_and_count;
    op_ptr->doWorkRhsHook = [&](DataOperator *base_op_ptr, int side,
                                EntityType type,
                                EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
      auto op_ptr = static_cast<VolOp *>(base_op_ptr);
      volume_and_count[VOL] += op_ptr->getMeasure();
      volume_and_count[COUNT] += 1;
      // in necessary at integration over Gauss points.
      MoFEMFunctionReturn(0);
    };
    volume_and_count = {0, 0};
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getOpPtrVector().push_back(op_ptr);
 
    auto dm = simple->getDM();
    CHK_THROW_MESSAGE(
        DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), fe),
        "cac volume");
    std::array<double, 2> tot_volume_and_count;
    MPI_Allreduce(volume_and_count.data(), tot_volume_and_count.data(),
                  volume_and_count.size(), MPI_DOUBLE, MPI_SUM,
                  mField.get_comm());
    return tot_volume_and_count;
  };
 
  meshVolumeAndCount = get_volume();
  MOFEM_LOG("PLASTICITY", Sev::inform)
      << "Mesh volume " << meshVolumeAndCount[VOL] << " nb. of elements "
      << meshVolumeAndCount[COUNT];
 
  MoFEMFunctionReturn(0);
}

setupProblem() [2/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [3/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [4/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [5/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [6/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [7/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [8/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [9/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [10/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [11/12]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [12/12]

MoFEMErrorCode Example::setupProblem ( )

private

Setup fields, approximation spaces and DOFs.

solveElastic()

MoFEMErrorCode Example::solveElastic ( )

private

Solve forward elastic problem.

[Push operators to pipeline]

[Solve]

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 1219 of file adjoint.cpp.

                                     {
  MoFEMFunctionBegin;
  auto simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
  auto f = createDMVector(dm);
  auto d = vectorDuplicate(f);
  CHKERR VecZeroEntries(d);
  CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
 
  auto set_essential_bc = [&]() {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
 
    auto ksp_ctx_ptr = getDMKspCtx(dm);
    auto pre_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_rhs = boost::make_shared<FEMethod>();
    auto post_proc_lhs = boost::make_shared<FEMethod>();
 
    auto get_pre_proc_hook = [&]() {
      return EssentialPreProc<DisplacementCubitBcData>(mField, pre_proc_rhs,
                                                       {});
    };
    pre_proc_rhs->preProcessHook = get_pre_proc_hook();
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(mField,
                                                           post_proc_rhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    auto get_post_proc_hook_lhs = [this, post_proc_lhs]() {
      MoFEMFunctionBeginHot;
 
      CHKERR EssentialPostProcLhs<DisplacementCubitBcData>(mField,
                                                           post_proc_lhs, 1.)();
      MoFEMFunctionReturnHot(0);
    };
 
    post_proc_rhs->postProcessHook = get_post_proc_hook_rhs;
    post_proc_lhs->postProcessHook = get_post_proc_hook_lhs;
 
    ksp_ctx_ptr->getPreProcComputeRhs().push_front(pre_proc_rhs);
    ksp_ctx_ptr->getPostProcComputeRhs().push_back(post_proc_rhs);
    ksp_ctx_ptr->getPostProcSetOperators().push_back(post_proc_lhs);
    MoFEMFunctionReturn(0);
  };
 
  auto setup_and_solve = [&](auto solver) {
    MoFEMFunctionBegin;
    BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp";
    CHKERR KSPSetUp(solver);
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSetUp <= Done";
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve";
    CHKERR KSPSolve(solver, f, d);
    MOFEM_LOG("TIMER", Sev::noisy) << "KSPSolve <= Done";
    MoFEMFunctionReturn(0);
  };
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
 
  CHKERR set_essential_bc();
  CHKERR setup_and_solve(kspElastic);
 
  CHKERR VecGhostUpdateBegin(d, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR VecGhostUpdateEnd(d, INSERT_VALUES, SCATTER_FORWARD);
  CHKERR DMoFEMMeshToLocalVector(dm, d, INSERT_VALUES, SCATTER_REVERSE);
 
  auto evaluate_field_at_the_point = [&]() {
    MoFEMFunctionBegin;
 
    int coords_dim = 3;
    std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
    PetscBool do_eval_field = PETSC_FALSE;
    CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                    field_eval_coords.data(), &coords_dim,
                                    &do_eval_field);
 
    if (do_eval_field) {
 
      vectorFieldPtr = boost::make_shared<MatrixDouble>();
      auto field_eval_data =
          mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->buildTree<SPACE_DIM>(field_eval_data, simple->getDomainFEName());
 
      field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
      auto no_rule = [](int, int, int) { return -1; };
      auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
      field_eval_fe_ptr->getRuleHook = no_rule;
 
      field_eval_fe_ptr->getOpPtrVector().push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("U", vectorFieldPtr));
 
      CHKERR mField.getInterface<FieldEvaluatorInterface>()
          ->evalFEAtThePoint<SPACE_DIM>(
              field_eval_coords.data(), 1e-12, simple->getProblemName(),
              simple->getDomainFEName(), field_eval_data,
              mField.get_comm_rank(), mField.get_comm_rank(), nullptr, MF_EXIST,
              QUIET);
 
      if (vectorFieldPtr->size1()) {
        auto t_disp = getFTensor1FromMat<SPACE_DIM>(*vectorFieldPtr);
        if constexpr (SPACE_DIM == 2)
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1);
        else
          MOFEM_LOG("FieldEvaluator", Sev::inform)
              << "U_X: " << t_disp(0) << " U_Y: " << t_disp(1)
              << " U_Z: " << t_disp(2);
      }
 
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR evaluate_field_at_the_point();
 
  MoFEMFunctionReturn(0);
}

solveSystem() [1/9]

MoFEMErrorCode Example::solveSystem ( )

private

[Solve]

[Push operators to pipeline]

[Solve]

< Mass matrix

< Linear solver

Examples

mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/mix-1/phase.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, and mofem/tutorials/vec-4/shallow_wave.cpp.

Definition at line 877 of file dynamic_first_order_con_law.cpp.

                                    {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  auto dm = simple->getDM();
 
  auto calculate_stress_ops = [&](auto &pip) {
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1});
 
    auto v_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("V", v_ptr));
    auto X_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", X_ptr));
 
    auto x_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("x_1", x_ptr));
 
    // Calculate unknown F
    auto mat_H_tensor_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateTensor2FieldValues<SPACE_DIM, SPACE_DIM>(
        "F", mat_H_tensor_ptr));
 
    auto u_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateDisplacement<SPACE_DIM>(x_ptr, X_ptr, u_ptr));
    // Calculate P
 
    auto mat_F_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateDeformationGradient<SPACE_DIM, SPACE_DIM>(
        mat_F_ptr, mat_H_tensor_ptr));
 
    PetscBool is_linear_elasticity = PETSC_TRUE;
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-is_linear_elasticity",
                               &is_linear_elasticity, PETSC_NULLPTR);
 
    auto mat_P_ptr = boost::make_shared<MatrixDouble>();
    if (is_linear_elasticity) {
      pip.push_back(new OpCalculatePiola<SPACE_DIM, SPACE_DIM>(
          shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_ptr,
          mat_F_ptr));
    } else {
      auto inv_F = boost::make_shared<MatrixDouble>();
      auto det_ptr = boost::make_shared<VectorDouble>();
 
      pip.push_back(new OpInvertMatrix<SPACE_DIM>(mat_F_ptr, det_ptr, inv_F));
 
      pip.push_back(new OpCalculatePiolaIncompressibleNH<SPACE_DIM, SPACE_DIM>(
          shear_modulus_G, bulk_modulus_K, mu, lamme_lambda, mat_P_ptr,
          mat_F_ptr, inv_F, det_ptr));
    }
 
    auto mat_v_grad_ptr = boost::make_shared<MatrixDouble>();
    pip.push_back(new OpCalculateVectorFieldGradient<SPACE_DIM, SPACE_DIM>(
        "V", mat_v_grad_ptr));
 
    return boost::make_tuple(v_ptr, X_ptr, x_ptr, mat_P_ptr, mat_F_ptr, u_ptr);
  };
 
  auto post_proc_boundary = [&]() {
    auto boundary_post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
 
    AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(
        boundary_post_proc_fe->getOpPtrVector(), {}, "GEOMETRY");
    auto op_loop_side =
        new OpLoopSide<SideEle>(mField, simple->getDomainFEName(), SPACE_DIM);
    // push ops to side element, through op_loop_side operator
    auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
          boundary_mat_F_ptr, boundary_u_ptr] =
        calculate_stress_ops(op_loop_side->getOpPtrVector());
    boundary_post_proc_fe->getOpPtrVector().push_back(op_loop_side);
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    boundary_post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            boundary_post_proc_fe->getPostProcMesh(),
            boundary_post_proc_fe->getMapGaussPts(),
 
            OpPPMap::DataMapVec{},
 
            OpPPMap::DataMapMat{{"V", boundary_v_ptr},
                                {"GEOMETRY", boundary_X_ptr},
                                {"x", boundary_x_ptr},
                                {"U", boundary_u_ptr}},
 
            OpPPMap::DataMapMat{{"FIRST_PIOLA", boundary_mat_P_ptr},
                                {"F", boundary_mat_F_ptr}},
 
            OpPPMap::DataMapMat{}
 
            )
 
    );
    return boundary_post_proc_fe;
  };
 
  // Add monitor to time solver
 
  double rho = 1.;
  CHKERR PetscOptionsGetReal(PETSC_NULLPTR, "", "-density", &rho, PETSC_NULLPTR);
  auto get_rho = [rho](const double, const double, const double) {
    return rho;
  };
 
  SmartPetscObj<Mat> M;   ///< Mass matrix
  SmartPetscObj<KSP> ksp; ///< Linear solver
 
  auto ts_pre_post_proc = boost::make_shared<TSPrePostProc>();
  tsPrePostProc = ts_pre_post_proc;
 
  CHKERR DMCreateMatrix_MoFEM(dm, M);
  CHKERR MatZeroEntries(M);
 
  boost::shared_ptr<DomainEle> vol_mass_ele(new DomainEle(mField));
 
  vol_mass_ele->B = M;
 
  auto integration_rule = [](int, int, int approx_order) {
    return 2 * approx_order;
  };
 
  vol_mass_ele->getRuleHook = integration_rule;
 
  vol_mass_ele->getOpPtrVector().push_back(new OpMassV("V", "V", get_rho));
  vol_mass_ele->getOpPtrVector().push_back(new OpMassF("F", "F"));
 
  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(), vol_mass_ele);
  CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
  CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 
  auto lumpVec = createDMVector(simple->getDM());
  CHKERR MatGetRowSum(M, lumpVec);
 
  CHKERR MatZeroEntries(M);
  CHKERR MatDiagonalSet(M, lumpVec, INSERT_VALUES);
 
  // Create and septup KSP (linear solver), we need this to calculate g(t,u) =
  // M^-1G(t,u)
  ksp = createKSP(mField.get_comm());
  CHKERR KSPSetOperators(ksp, M, M);
  CHKERR KSPSetFromOptions(ksp);
  CHKERR KSPSetUp(ksp);
 
  auto solve_boundary_for_g = [&]() {
    MoFEMFunctionBegin;
    if (*(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch)) {
 
      CHKERR VecGhostUpdateBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
                                 ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F,
                               ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecAssemblyBegin(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
      CHKERR VecAssemblyEnd(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
      *(pipeline_mng->getBoundaryExplicitRhsFE()->vecAssembleSwitch) = false;
 
      auto D =
          smartVectorDuplicate(pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
      CHKERR KSPSolve(ksp, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F, D);
      CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecCopy(D, pipeline_mng->getBoundaryExplicitRhsFE()->ts_F);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  pipeline_mng->getBoundaryExplicitRhsFE()->postProcessHook =
      solve_boundary_for_g;
 
  MoFEM::SmartPetscObj<TS> ts;
  ts = pipeline_mng->createTSEX(dm);
 
  // Field eval
  PetscBool field_eval_flag = PETSC_TRUE;
  boost::shared_ptr<MatrixDouble> velocity_field_ptr;
  boost::shared_ptr<MatrixDouble> geometry_field_ptr;
  boost::shared_ptr<MatrixDouble> spatial_position_field_ptr;
  boost::shared_ptr<SetPtsData> field_eval_data;
 
  std::array<double, 3> field_eval_coords = {0.5, 0.5, 5.};
  int dim = 3;
  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                  field_eval_coords.data(), &dim,
                                  &field_eval_flag);
 
  if (field_eval_flag) {
    field_eval_data =
        mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
    CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
        field_eval_data, simple->getDomainFEName());
 
    field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
 
    auto no_rule = [](int, int, int) { return -1; };
 
    auto fe_ptr = field_eval_data->feMethodPtr.lock();
    fe_ptr->getRuleHook = no_rule;
    velocity_field_ptr = boost::make_shared<MatrixDouble>();
    geometry_field_ptr = boost::make_shared<MatrixDouble>();
    spatial_position_field_ptr = boost::make_shared<MatrixDouble>();
    fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("V", velocity_field_ptr));
    fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY",
                                                    geometry_field_ptr));
    fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>(
            "x_2", spatial_position_field_ptr));
  }
 
  auto post_proc_domain = [&]() {
    auto post_proc_fe_vol = boost::make_shared<PostProcEle>(mField);
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    auto [boundary_v_ptr, boundary_X_ptr, boundary_x_ptr, boundary_mat_P_ptr,
          boundary_mat_F_ptr, boundary_u_ptr] =
        calculate_stress_ops(post_proc_fe_vol->getOpPtrVector());
 
    post_proc_fe_vol->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe_vol->getPostProcMesh(),
            post_proc_fe_vol->getMapGaussPts(),
 
            {},
 
            {{"V", boundary_v_ptr},
             {"GEOMETRY", boundary_X_ptr},
             {"x", boundary_x_ptr},
             {"U", boundary_u_ptr}},
 
            {{"FIRST_PIOLA", boundary_mat_P_ptr}, {"F", boundary_mat_F_ptr}},
 
            {}
 
            )
 
    );
    return post_proc_fe_vol;
  };
 
  boost::shared_ptr<FEMethod> null_fe;
  auto monitor_ptr = boost::make_shared<Monitor>(
      SmartPetscObj<DM>(dm, true), mField, post_proc_domain(),
      post_proc_boundary(), velocity_field_ptr, spatial_position_field_ptr,
      geometry_field_ptr, field_eval_coords, field_eval_data);
 
  CHKERR DMMoFEMTSSetMonitor(dm, ts, simple->getDomainFEName(), null_fe,
                             null_fe, monitor_ptr);
 
  double ftime = 1;
  // CHKERR TSSetMaxTime(ts, ftime);
  CHKERR TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
 
  auto T = createDMVector(simple->getDM());
  CHKERR DMoFEMMeshToLocalVector(simple->getDM(), T, INSERT_VALUES,
                                 SCATTER_FORWARD);
  CHKERR TSSetSolution(ts, T);
  CHKERR TSSetFromOptions(ts);
 
  CHKERR TSSetPostStage(ts, TSPrePostProc::tsPostStage);
  CHKERR TSSetPostStep(ts, TSPrePostProc::tsPostStep);
  CHKERR TSSetPreStep(ts, TSPrePostProc::tsPreStep);
 
  boost::shared_ptr<ForcesAndSourcesCore> null;
 
  if (auto ptr = tsPrePostProc.lock()) {
    ptr->fsRawPtr = this;
    CHKERR TSSetUp(ts);
    CHKERR TSSolve(ts, NULL);
    CHKERR TSGetTime(ts, &ftime);
  }
 
  MoFEMFunctionReturn(0);
}

solveSystem() [2/9]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [3/9]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [4/9]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [5/9]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [6/9]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [7/9]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [8/9]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [9/9]

MoFEMErrorCode Example::solveSystem ( )

private

testOperators()

MoFEMErrorCode Example::testOperators ( )

private

[Solve]

[TestOperators]

Examples

mofem/tutorials/adv-0/plastic.cpp.

Definition at line 1477 of file plastic.cpp.

                                      {
  MoFEMFunctionBegin;
 
  // get operators tester
  auto simple = mField.getInterface<Simple>();
  auto opt = mField.getInterface<OperatorsTester>(); // get interface to
                                                     // OperatorsTester
  auto pip = mField.getInterface<PipelineManager>(); // get interface to
                                                     // pipeline manager
 
  constexpr double eps = 1e-9;
 
  auto x = opt->setRandomFields(simple->getDM(), {
 
                                                     {"U", {-1e-6, 1e-6}},
 
                                                     {"EP", {-1e-6, 1e-6}},
 
                                                     {"TAU", {0, 1e-4}}
 
                                                 });
 
  auto dot_x_plastic_active =
      opt->setRandomFields(simple->getDM(), {
 
                                                {"U", {-1, 1}},
 
                                                {"EP", {-1, 1}},
 
                                                {"TAU", {0.1, 1}}
 
                                            });
  auto diff_x_plastic_active =
      opt->setRandomFields(simple->getDM(), {
 
                                                {"U", {-1, 1}},
 
                                                {"EP", {-1, 1}},
 
                                                {"TAU", {-1, 1}}
 
                                            });
 
  auto dot_x_elastic =
      opt->setRandomFields(simple->getDM(), {
 
                                                {"U", {-1, 1}},
 
                                                {"EP", {-1, 1}},
 
                                                {"TAU", {-1, -0.1}}
 
                                            });
  auto diff_x_elastic =
      opt->setRandomFields(simple->getDM(), {
 
                                                {"U", {-1, 1}},
 
                                                {"EP", {-1, 1}},
 
                                                {"TAU", {-1, 1}}
 
                                            });
 
  auto test_domain_ops = [&](auto fe_name, auto lhs_pipeline, auto rhs_pipeline,
                             auto dot_x, auto diff_x) {
    MoFEMFunctionBegin;
 
    auto diff_res = opt->checkCentralFiniteDifference(
        simple->getDM(), fe_name, rhs_pipeline, lhs_pipeline, x, dot_x,
        SmartPetscObj<Vec>(), diff_x, 0, 0.5, eps);
 
    // Calculate norm of difference between directional derivative calculated
    // from finite difference, and tangent matrix.
    double fnorm;
    CHKERR VecNorm(diff_res, NORM_2, &fnorm);
    MOFEM_LOG_C("PLASTICITY", Sev::inform,
                "Test consistency of tangent matrix %3.4e", fnorm);
 
    constexpr double err = 1e-5;
    if (fnorm > err)
      SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
               "Norm of directional derivative too large err = %3.4e", fnorm);
 
    MoFEMFunctionReturn(0);
  };
 
  MOFEM_LOG("PLASTICITY", Sev::inform) << "Elastic active";
  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
                         pip->getDomainRhsFE(), dot_x_elastic, diff_x_elastic);
 
  MOFEM_LOG("PLASTICITY", Sev::inform) << "Plastic active";
  CHKERR test_domain_ops(simple->getDomainFEName(), pip->getDomainLhsFE(),
                         pip->getDomainRhsFE(), dot_x_plastic_active,
                         diff_x_plastic_active);
 
  MoFEMFunctionReturn(0);
};

topologyModes()

MoFEMErrorCode Example::topologyModes ( )

private

Compute topology optimization modes.

[Boundary condition]

[Adjoint modes]

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 746 of file adjoint.cpp.

                                      {
  MoFEMFunctionBegin;
 
  auto opt_ents = get_range_from_block(mField, "OPTIMISE", SPACE_DIM - 1);
  auto subset_dm_bdy = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMSetSquareProblem(subset_dm_bdy, PETSC_TRUE);
  CHKERR DMMoFEMCreateSubDM(subset_dm_bdy, adjointDM, "SUBSET_BDY");
  CHKERR DMMoFEMAddElement(subset_dm_bdy, "ADJOINT_BOUNDARY_FE");
  CHKERR DMMoFEMAddSubFieldRow(subset_dm_bdy, "ADJOINT_FIELD",
                               boost::make_shared<Range>(opt_ents));
  CHKERR DMMoFEMAddSubFieldCol(subset_dm_bdy, "ADJOINT_FIELD",
                               boost::make_shared<Range>(opt_ents));
  CHKERR DMSetUp(subset_dm_bdy);
 
  auto subset_dm_domain = createDM(mField.get_comm(), "DMMOFEM");
  CHKERR DMMoFEMSetSquareProblem(subset_dm_domain, PETSC_TRUE);
  CHKERR DMMoFEMCreateSubDM(subset_dm_domain, adjointDM, "SUBSET_DOMAIN");
  CHKERR DMMoFEMAddElement(subset_dm_domain, "ADJOINT_DOMAIN_FE");
  CHKERR DMMoFEMAddSubFieldRow(subset_dm_domain, "ADJOINT_FIELD");
  CHKERR DMMoFEMAddSubFieldCol(subset_dm_domain, "ADJOINT_FIELD");
  CHKERR DMSetUp(subset_dm_domain);
 
  // remove dofs on boundary of the domain
  auto remove_dofs = [&]() {
    MoFEMFunctionBegin;
 
    std::array<Range, 3> remove_dim_ents;
    remove_dim_ents[0] =
        get_range_from_block(mField, "OPT_REMOVE_X", SPACE_DIM - 1);
    remove_dim_ents[1] =
        get_range_from_block(mField, "OPT_REMOVE_Y", SPACE_DIM - 1);
    remove_dim_ents[2] =
        get_range_from_block(mField, "OPT_REMOVE_Z", SPACE_DIM - 1);
 
    for (int d = 0; d != 3; ++d) {
      MOFEM_LOG("WORLD", Sev::inform)
          << "Removing topology modes on block OPT_REMOVE_" << (char)('X' + d)
          << " with " << remove_dim_ents[d].size() << " entities";
    }
 
    Range body_ents;
    CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, body_ents,
                                                       true);
    auto skin = moab::Skinner(&mField.get_moab());
    Range boundary_ents;
    CHKERR skin.find_skin(0, body_ents, false, boundary_ents);
    for (int d = 0; d != 3; ++d) {
      boundary_ents = subtract(boundary_ents, remove_dim_ents[d]);
    }
    ParallelComm *pcomm =
        ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
    CHKERR pcomm->filter_pstatus(boundary_ents,
                                 PSTATUS_SHARED | PSTATUS_MULTISHARED,
                                 PSTATUS_NOT, -1, nullptr);
    for (auto d = SPACE_DIM - 2; d >= 0; --d) {
      if (d >= 0) {
        Range ents;
        CHKERR mField.get_moab().get_adjacencies(boundary_ents, d, false, ents,
                                                 moab::Interface::UNION);
        boundary_ents.merge(ents);
      } else {
        Range verts;
        CHKERR mField.get_moab().get_connectivity(boundary_ents, verts);
        boundary_ents.merge(verts);
      }
      CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
          boundary_ents);
    }
    boundary_ents.merge(opt_ents);
    CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
        "SUBSET_DOMAIN", "ADJOINT_FIELD", boundary_ents);
    for (int d = 0; d != 3; ++d) {
      CHKERR mField.getInterface<ProblemsManager>()->removeDofsOnEntities(
          "SUBSET_DOMAIN", "ADJOINT_FIELD", remove_dim_ents[d], d, d);
    }
 
    // #ifndef NDEBUG
    if (mField.get_comm_rank() == 0) {
      CHKERR save_range(mField.get_moab(), "topoMode_boundary_ents.vtk",
                        boundary_ents);
    }
    // #endif
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR remove_dofs();
 
  auto get_lhs_fe = [&]() {
    auto fe_lhs = boost::make_shared<BoundaryEle>(mField);
    fe_lhs->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_lhs->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
                                                              "GEOMETRY");
    using OpMass = FormsIntegrators<BoundaryEleOp>::Assembly<A>::BiLinearForm<
        I>::OpMass<1, SPACE_DIM>;
    pip.push_back(new OpMass("ADJOINT_FIELD", "ADJOINT_FIELD",
                             [](double, double, double) { return 1.; }));
    return fe_lhs;
  };
 
  auto get_rhs_fe = [&]() {
    auto fe_rhs = boost::make_shared<BoundaryEle>(mField);
    fe_rhs->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe_rhs->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM - 1, SPACE_DIM, SPACE_DIM>::add(pip, {},
                                                              "GEOMETRY");
 
    return fe_rhs;
  };
 
  auto block_name = "OPTIMISE";
  auto mesh_mng = mField.getInterface<MeshsetsManager>();
  auto bcs = mesh_mng->getCubitMeshsetPtr(
 
      std::regex((boost::format("%s(.*)") % block_name).str())
 
  );
 
  for (auto &v : modeVecs) {
    v = createDMVector(subset_dm_bdy);
  }
 
  using OP = FormsIntegrators<BoundaryEleOp>::Assembly<A>::OpBase;
  struct OpMode : public OP {
    OpMode(const std::string name,
           boost::shared_ptr<ObjectiveFunctionData> python_ptr, int id,
           std::vector<SmartPetscObj<Vec>> mode_vecs,
           std::vector<std::array<double, 3>> mode_centroids,
           std::vector<std::array<double, 6>> mode_bboxes, int block_counter,
           int mode_counter, boost::shared_ptr<Range> range = nullptr)
        : OP(name, name, OP::OPROW, range), pythonPtr(python_ptr), iD(id),
          modeVecs(mode_vecs), modeCentroids(mode_centroids),
          modeBboxes(mode_bboxes), blockCounter(block_counter),
          modeCounter(mode_counter) {}
 
    MoFEMErrorCode doWork(int side, EntityType type, EntData &data) {
      MoFEMFunctionBegin;
 
      if (OP::entsPtr) {
        if (OP::entsPtr->find(this->getFEEntityHandle()) == OP::entsPtr->end())
          MoFEMFunctionReturnHot(0);
      }
 
      auto nb_rows = data.getIndices().size();
      if (!nb_rows) {
        MoFEMFunctionReturnHot(0);
      }
      auto nb_base_functions = data.getN().size2();
 
      FTENSOR_INDEX(SPACE_DIM, i);
      CHKERR pythonPtr->blockModes(iD, OP::getCoordsAtGaussPts(),
                                   modeCentroids[blockCounter],
                                   modeBboxes[blockCounter], blockModes);
 
      auto nb_integration_pts = getGaussPts().size2();
      if (blockModes.size2() != 3 * nb_integration_pts) {
        MOFEM_LOG("WORLD", Sev::error)
            << "Number of modes does not match number of integration points: "
            << blockModes.size2() << "!=" << 3 * nb_integration_pts;
        CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "modes/integration points");
      }
 
      VectorDouble nf(nb_rows);
 
      int nb_modes = blockModes.size1();
      for (auto mode = 0; mode != nb_modes; ++mode) {
        nf.clear();
        // get mode
        auto t_mode = getFTensor1FromPtr<3>(&blockModes(mode, 0));
        // get element volume
        const double vol = OP::getMeasure();
        // get integration weights
        auto t_w = OP::getFTensor0IntegrationWeight();
        // get base function gradient on rows
        auto t_base = data.getFTensor0N();
        // loop over integration points
        for (int gg = 0; gg != nb_integration_pts; gg++) {
 
          // take into account Jacobian
          const double alpha = t_w * vol;
          // loop over rows base functions
          auto t_nf = getFTensor1FromPtr<SPACE_DIM>(nf.data().data());
          int rr = 0;
          for (; rr != nb_rows / SPACE_DIM; ++rr) {
            t_nf(i) += alpha * t_base * t_mode(i);
            ++t_base;
            ++t_nf;
          }
          for (; rr < nb_base_functions; ++rr)
            ++t_base;
          ++t_w;    // move to another integration weight
          ++t_mode; // move to another mode
        }
        Vec vec = modeVecs[modeCounter + mode];
        auto size = data.getIndices().size();
        auto *indices = data.getIndices().data().data();
        auto *nf_data = nf.data().data();
        CHKERR VecSetValues(vec, size, indices, nf_data, ADD_VALUES);
      }
 
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<ObjectiveFunctionData> pythonPtr;
    MatrixDouble blockModes;
    std::vector<std::array<double, 3>> modeCentroids;
    std::vector<std::array<double, 6>> modeBboxes;
    int iD;
    std::vector<SmartPetscObj<Vec>> modeVecs;
    int blockCounter;
    int modeCounter;
  };
 
  auto solve_bdy = [&]() {
    MoFEMFunctionBegin;
 
    auto fe_lhs = get_lhs_fe();
    auto fe_rhs = get_rhs_fe();
    int block_counter = 0;
    int mode_counter = 0;
    for (auto &bc : bcs) {
      auto id = bc->getMeshsetId();
      Range ents;
      CHKERR mField.get_moab().get_entities_by_handle(bc->getMeshset(), ents,
                                                      true);
      auto range = boost::make_shared<Range>(ents);
      auto &pip_rhs = fe_rhs->getOpPtrVector();
      pip_rhs.push_back(new OpMode("ADJOINT_FIELD", pythonPtr, id, modeVecs,
                                   modeCentroids, modeBBoxes, block_counter,
                                   mode_counter, range));
      CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
                                      fe_rhs);
      pip_rhs.pop_back();
      int nb_modes;
      CHKERR pythonPtr->numberOfModes(id, nb_modes);
      ++block_counter;
      mode_counter += nb_modes;
      MOFEM_LOG("WORLD", Sev::inform)
          << "Setting mode block block: " << bc->getName()
          << " with ID: " << bc->getMeshsetId()
          << " total modes: " << mode_counter;
    }
 
    for (auto &v : modeVecs) {
      CHKERR VecAssemblyBegin(v);
      CHKERR VecAssemblyEnd(v);
      CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
    }
 
    auto M = createDMMatrix(subset_dm_bdy);
    fe_lhs->B = M;
    CHKERR DMoFEMLoopFiniteElements(subset_dm_bdy, "ADJOINT_BOUNDARY_FE",
                                    fe_lhs);
    CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 
    auto solver = createKSP(mField.get_comm());
    CHKERR KSPSetOperators(solver, M, M);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
    auto v = createDMVector(subset_dm_bdy);
    for (auto &f : modeVecs) {
      CHKERR KSPSolve(solver, f, v);
      CHKERR VecSwap(f, v);
    }
 
    for (auto &v : modeVecs) {
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
    }
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_bdy();
 
  auto get_elastic_fe_lhs = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    CHKERR HookeOps::opFactoryDomainLhs<SPACE_DIM, A, I, DomainEleOp>(
        mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
    return fe;
  };
 
  auto get_elastic_fe_rhs = [&]() {
    auto fe = boost::make_shared<DomainEle>(mField);
    fe->getRuleHook = [](int, int, int p_data) {
      return 2 * p_data + p_data - 1;
    };
    auto &pip = fe->getOpPtrVector();
    CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1},
                                                          "GEOMETRY");
    CHKERR HookeOps::opFactoryDomainRhs<SPACE_DIM, A, I, DomainEleOp>(
        mField, pip, "ADJOINT_FIELD", "MAT_ADJOINT", Sev::noisy);
    return fe;
  };
 
  auto adjoint_gradient_postprocess = [&](auto mode) {
    MoFEMFunctionBegin;
    auto post_proc_mesh = boost::make_shared<moab::Core>();
    auto post_proc_begin =
        boost::make_shared<PostProcBrokenMeshInMoabBaseBegin>(mField,
                                                              post_proc_mesh);
    auto post_proc_end = boost::make_shared<PostProcBrokenMeshInMoabBaseEnd>(
        mField, post_proc_mesh);
 
    auto geom_vec = boost::make_shared<MatrixDouble>();
 
    auto post_proc_fe =
        boost::make_shared<PostProcEleDomain>(mField, post_proc_mesh);
    AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        post_proc_fe->getOpPtrVector(), {H1}, "GEOMETRY");
    post_proc_fe->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("ADJOINT_FIELD", geom_vec,
                                                    modeVecs[mode]));
 
    using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
    post_proc_fe->getOpPtrVector().push_back(
 
        new OpPPMap(
 
            post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
 
            {},
 
            {{"MODE", geom_vec}},
 
            {},
 
            {}
 
            )
 
    );
 
    CHKERR DMoFEMPreProcessFiniteElements(adjointDM,
                                          post_proc_begin->getFEMethod());
    CHKERR DMoFEMLoopFiniteElements(adjointDM, "ADJOINT_DOMAIN_FE",
                                    post_proc_fe);
    CHKERR DMoFEMPostProcessFiniteElements(adjointDM,
                                           post_proc_begin->getFEMethod());
 
    CHKERR post_proc_end->writeFile("mode_" + std::to_string(mode) + ".h5m");
 
    MoFEMFunctionReturn(0);
  };
 
  auto solve_domain = [&]() {
    MoFEMFunctionBegin;
    auto fe_lhs = get_elastic_fe_lhs();
    auto fe_rhs = get_elastic_fe_rhs();
    auto v = createDMVector(subset_dm_domain);
    auto F = vectorDuplicate(v);
    fe_rhs->f = F;
 
    auto M = createDMMatrix(subset_dm_domain);
    fe_lhs->B = M;
    CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
                                    fe_lhs);
    CHKERR MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);
    CHKERR MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);
 
    auto solver = createKSP(mField.get_comm());
    CHKERR KSPSetOperators(solver, M, M);
    CHKERR KSPSetFromOptions(solver);
    CHKERR KSPSetUp(solver);
 
    int mode_counter = 0;
    for (auto &f : modeVecs) {
      CHKERR mField.getInterface<FieldBlas>()->setField(0, "ADJOINT_FIELD");
      CHKERR DMoFEMMeshToLocalVector(subset_dm_bdy, f, INSERT_VALUES,
                                     SCATTER_REVERSE);
      CHKERR VecZeroEntries(F);
      CHKERR DMoFEMLoopFiniteElements(subset_dm_domain, "ADJOINT_DOMAIN_FE",
                                      fe_rhs);
      CHKERR VecAssemblyBegin(F);
      CHKERR VecAssemblyEnd(F);
      CHKERR VecGhostUpdateBegin(F, ADD_VALUES, SCATTER_REVERSE);
      CHKERR VecGhostUpdateEnd(F, ADD_VALUES, SCATTER_REVERSE);
      CHKERR KSPSolve(solver, F, v);
      CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
      CHKERR DMoFEMMeshToLocalVector(subset_dm_domain, v, INSERT_VALUES,
                                     SCATTER_REVERSE);
      auto m = createDMVector(adjointDM);
      CHKERR DMoFEMMeshToLocalVector(adjointDM, m, INSERT_VALUES,
                                     SCATTER_FORWARD);
      f = m;
      ++mode_counter;
    }
    MoFEMFunctionReturn(0);
  };
 
  CHKERR solve_domain();
 
  for (int i = 0; i < modeVecs.size(); ++i) {
    CHKERR adjoint_gradient_postprocess(i);
  }
 
  MoFEMFunctionReturn(0);
}

tsSolve()

MoFEMErrorCode Example::tsSolve ( )

private

Examples

mofem/tutorials/adv-0/plastic.cpp.

Definition at line 834 of file plastic.cpp.

                                {
  MoFEMFunctionBegin;
 
  Simple *simple = mField.getInterface<Simple>();
  PipelineManager *pip_mng = mField.getInterface<PipelineManager>();
  ISManager *is_manager = mField.getInterface<ISManager>();
 
  auto snes_ctx_ptr = getDMSnesCtx(simple->getDM());
 
  auto set_section_monitor = [&](auto solver) {
    MoFEMFunctionBegin;
    SNES snes;
    CHKERR TSGetSNES(solver, &snes);
    CHKERR SNESMonitorSet(snes,
                          (MoFEMErrorCode(*)(SNES, PetscInt, PetscReal,
                                             void *))MoFEMSNESMonitorFields,
                          (void *)(snes_ctx_ptr.get()), nullptr);
    MoFEMFunctionReturn(0);
  };
 
  auto create_post_process_elements = [&]() {
    auto push_vol_ops = [this](auto &pip) {
      CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
          pip, {H1, HDIV}, "GEOMETRY");
 
      auto [common_plastic_ptr, common_hencky_ptr] =
          PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
              mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
 
      if (common_hencky_ptr) {
        if (common_plastic_ptr->mGradPtr != common_hencky_ptr->matGradPtr)
          CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
      }
 
      return std::make_pair(common_plastic_ptr, common_hencky_ptr);
    };
 
    auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
      MoFEMFunctionBegin;
 
      auto &pip = pp_fe->getOpPtrVector();
 
      auto [common_plastic_ptr, common_hencky_ptr] = p;
 
      using OpPPMap = OpPostProcMapInMoab<SPACE_DIM, SPACE_DIM>;
 
      auto x_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(
          new OpCalculateVectorFieldValues<SPACE_DIM>("GEOMETRY", x_ptr));
      auto u_ptr = boost::make_shared<MatrixDouble>();
      pip.push_back(new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
      if (is_large_strains) {
 
        pip.push_back(
 
            new OpPPMap(
 
                pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
                {{"PLASTIC_SURFACE",
                  common_plastic_ptr->getPlasticSurfacePtr()},
                 {"PLASTIC_MULTIPLIER",
                  common_plastic_ptr->getPlasticTauPtr()}},
 
                {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
 
                {{"GRAD", common_hencky_ptr->matGradPtr},
                 {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()}},
 
                {{"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()},
                 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
                 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
 
                )
 
        );
 
      } else {
 
        pip.push_back(
 
            new OpPPMap(
 
                pp_fe->getPostProcMesh(), pp_fe->getMapGaussPts(),
 
                {{"PLASTIC_SURFACE",
                  common_plastic_ptr->getPlasticSurfacePtr()},
                 {"PLASTIC_MULTIPLIER",
                  common_plastic_ptr->getPlasticTauPtr()}},
 
                {{"U", u_ptr}, {"GEOMETRY", x_ptr}},
 
                {},
 
                {{"STRAIN", common_plastic_ptr->mStrainPtr},
                 {"STRESS", common_plastic_ptr->mStressPtr},
                 {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()},
                 {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()}}
 
                )
 
        );
      }
 
      MoFEMFunctionReturn(0);
    };
 
    PetscBool post_proc_vol;
    PetscBool post_proc_skin;
 
    if constexpr (SPACE_DIM == 2) {
      post_proc_vol = PETSC_TRUE;
      post_proc_skin = PETSC_FALSE;
    } else {
      post_proc_vol = PETSC_FALSE;
      post_proc_skin = PETSC_TRUE;
    }
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_vol", &post_proc_vol,
                               PETSC_NULLPTR);
    CHKERR PetscOptionsGetBool(PETSC_NULLPTR, "", "-post_proc_skin",
                               &post_proc_skin, PETSC_NULLPTR);
 
    auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
                          post_proc_vol]() {
      if (post_proc_vol == PETSC_FALSE)
        return boost::shared_ptr<PostProcEle>();
      auto pp_fe = boost::make_shared<PostProcEle>(mField);
      CHK_MOAB_THROW(
          push_vol_post_proc_ops(pp_fe, push_vol_ops(pp_fe->getOpPtrVector())),
          "push_vol_post_proc_ops");
      return pp_fe;
    };
 
    auto skin_post_proc = [this, push_vol_post_proc_ops, push_vol_ops,
                           post_proc_skin]() {
      if (post_proc_skin == PETSC_FALSE)
        return boost::shared_ptr<SkinPostProcEle>();
 
      auto simple = mField.getInterface<Simple>();
      auto pp_fe = boost::make_shared<SkinPostProcEle>(mField);
      auto op_side = new OpLoopSide<SideEle>(mField, simple->getDomainFEName(),
                                             SPACE_DIM, Sev::verbose);
      pp_fe->getOpPtrVector().push_back(op_side);
      CHK_MOAB_THROW(push_vol_post_proc_ops(
                         pp_fe, push_vol_ops(op_side->getOpPtrVector())),
                     "push_vol_post_proc_ops");
      return pp_fe;
    };
 
    return std::make_pair(vol_post_proc(), skin_post_proc());
  };
 
  auto scatter_create = [&](auto D, auto coeff) {
    SmartPetscObj<IS> is;
    CHKERR is_manager->isCreateProblemFieldAndRank(simple->getProblemName(),
                                                   ROW, "U", coeff, coeff, is);
    int loc_size;
    CHKERR ISGetLocalSize(is, &loc_size);
    Vec v;
    CHKERR VecCreateMPI(mField.get_comm(), loc_size, PETSC_DETERMINE, &v);
    VecScatter scatter;
    CHKERR VecScatterCreate(D, is, v, PETSC_NULLPTR, &scatter);
    return std::make_tuple(SmartPetscObj<Vec>(v),
                           SmartPetscObj<VecScatter>(scatter));
  };
 
  boost::shared_ptr<SetPtsData> field_eval_data;
  boost::shared_ptr<MatrixDouble> u_field_ptr;
 
  std::array<double, 3> field_eval_coords{0.0, 0.0, 0.0};
  int coords_dim = 3;
  CHKERR PetscOptionsGetRealArray(NULL, NULL, "-field_eval_coords",
                                  field_eval_coords.data(), &coords_dim,
                                  &do_eval_field);
 
  boost::shared_ptr<std::map<std::string, boost::shared_ptr<VectorDouble>>>
      scalar_field_ptrs = boost::make_shared<
          std::map<std::string, boost::shared_ptr<VectorDouble>>>();
  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
      vector_field_ptrs = boost::make_shared<
          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
      sym_tensor_field_ptrs = boost::make_shared<
          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
  boost::shared_ptr<std::map<std::string, boost::shared_ptr<MatrixDouble>>>
      tensor_field_ptrs = boost::make_shared<
          std::map<std::string, boost::shared_ptr<MatrixDouble>>>();
 
  if (do_eval_field) {
    auto u_field_ptr = boost::make_shared<MatrixDouble>();
    field_eval_data =
        mField.getInterface<FieldEvaluatorInterface>()->getData<DomainEle>();
 
    CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree<SPACE_DIM>(
        field_eval_data, simple->getDomainFEName());
 
    field_eval_data->setEvalPoints(field_eval_coords.data(), 1);
    auto no_rule = [](int, int, int) { return -1; };
    auto field_eval_fe_ptr = field_eval_data->feMethodPtr.lock();
    field_eval_fe_ptr->getRuleHook = no_rule;
 
    CHKERR PlasticOps::AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
        field_eval_fe_ptr->getOpPtrVector(), {H1, HDIV}, "GEOMETRY");
 
    auto [common_plastic_ptr, common_hencky_ptr] =
        PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
            mField, "MAT_PLASTIC", field_eval_fe_ptr->getOpPtrVector(), "U",
            "EP", "TAU", 1., Sev::inform);
 
    field_eval_fe_ptr->getOpPtrVector().push_back(
        new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_field_ptr));
 
    if ((common_plastic_ptr) && (common_hencky_ptr) && (scalar_field_ptrs)) {
      if (is_large_strains) {
        scalar_field_ptrs->insert(
            {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
        scalar_field_ptrs->insert(
            {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
        vector_field_ptrs->insert({"U", u_field_ptr});
        sym_tensor_field_ptrs->insert(
            {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
        sym_tensor_field_ptrs->insert(
            {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
        sym_tensor_field_ptrs->insert(
            {"HENCKY_STRAIN", common_hencky_ptr->getMatLogC()});
        tensor_field_ptrs->insert({"GRAD", common_hencky_ptr->matGradPtr});
        tensor_field_ptrs->insert(
            {"FIRST_PIOLA", common_hencky_ptr->getMatFirstPiolaStress()});
      } else {
        scalar_field_ptrs->insert(
            {"PLASTIC_SURFACE", common_plastic_ptr->getPlasticSurfacePtr()});
        scalar_field_ptrs->insert(
            {"PLASTIC_MULTIPLIER", common_plastic_ptr->getPlasticTauPtr()});
        vector_field_ptrs->insert({"U", u_field_ptr});
        sym_tensor_field_ptrs->insert(
            {"STRAIN", common_plastic_ptr->mStrainPtr});
        sym_tensor_field_ptrs->insert(
            {"STRESS", common_plastic_ptr->mStressPtr});
        sym_tensor_field_ptrs->insert(
            {"PLASTIC_STRAIN", common_plastic_ptr->getPlasticStrainPtr()});
        sym_tensor_field_ptrs->insert(
            {"PLASTIC_FLOW", common_plastic_ptr->getPlasticFlowPtr()});
      }
    }
  }
 
  auto test_monitor_ptr = boost::make_shared<FEMethod>();
 
  auto set_time_monitor = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    boost::shared_ptr<Monitor<SPACE_DIM>> monitor_ptr(new Monitor<SPACE_DIM>(
        dm, create_post_process_elements(), reactionFe, uXScatter, uYScatter,
        uZScatter, field_eval_coords, field_eval_data, scalar_field_ptrs,
        vector_field_ptrs, sym_tensor_field_ptrs, tensor_field_ptrs));
    boost::shared_ptr<ForcesAndSourcesCore> null;
 
    test_monitor_ptr->postProcessHook = [&]() {
      MoFEMFunctionBegin;
 
      if (atom_test && fabs(test_monitor_ptr->ts_t - 0.5) < 1e-12 &&
          test_monitor_ptr->ts_step == 25) {
 
        if (scalar_field_ptrs->at("PLASTIC_MULTIPLIER")->size()) {
          auto t_tau =
              getFTensor0FromVec(*scalar_field_ptrs->at("PLASTIC_MULTIPLIER"));
          MOFEM_LOG("PlasticSync", Sev::inform) << "Eval point tau: " << t_tau;
 
          if (atom_test == 1 && fabs(t_tau - 0.688861) > 1e-5) {
            SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                     "atom test %d failed: wrong plastic multiplier value",
                     atom_test);
          }
        }
 
        if (vector_field_ptrs->at("U")->size1()) {
          FTensor::Index<'i', SPACE_DIM> i;
          auto t_disp =
              getFTensor1FromMat<SPACE_DIM>(*vector_field_ptrs->at("U"));
          MOFEM_LOG("PlasticSync", Sev::inform) << "Eval point U: " << t_disp;
 
          if (atom_test == 1 && fabs(t_disp(0) - 0.25 / 2.) > 1e-5 ||
              fabs(t_disp(1) + 0.0526736) > 1e-5) {
            SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                     "atom test %d failed: wrong displacement value",
                     atom_test);
          }
        }
 
        if (sym_tensor_field_ptrs->at("PLASTIC_STRAIN")->size1()) {
          auto t_plastic_strain = getFTensor2SymmetricFromMat<SPACE_DIM>(
              *sym_tensor_field_ptrs->at("PLASTIC_STRAIN"));
          MOFEM_LOG("PlasticSync", Sev::inform)
              << "Eval point EP: " << t_plastic_strain;
 
          if (atom_test == 1 &&
                  fabs(t_plastic_strain(0, 0) - 0.221943) > 1e-5 ||
              fabs(t_plastic_strain(0, 1)) > 1e-5 ||
              fabs(t_plastic_strain(1, 1) + 0.110971) > 1e-5) {
            SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                     "atom test %d failed: wrong plastic strain value",
                     atom_test);
          }
        }
 
        if (tensor_field_ptrs->at("FIRST_PIOLA")->size1()) {
          auto t_piola_stress = getFTensor2FromMat<SPACE_DIM, SPACE_DIM>(
              *tensor_field_ptrs->at("FIRST_PIOLA"));
          MOFEM_LOG("PlasticSync", Sev::inform)
              << "Eval point Piola stress: " << t_piola_stress;
 
          if (atom_test == 1 && fabs((t_piola_stress(0, 0) - 198.775) /
                                     t_piola_stress(0, 0)) > 1e-5 ||
              fabs(t_piola_stress(0, 1)) + fabs(t_piola_stress(1, 0)) +
                      fabs(t_piola_stress(1, 1)) >
                  1e-5) {
            SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
                     "atom test %d failed: wrong Piola stress value",
                     atom_test);
          }
        }
      }
 
      MOFEM_LOG_SYNCHRONISE(mField.get_comm());
      MoFEMFunctionReturn(0);
    };
 
    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                               monitor_ptr, null, test_monitor_ptr);
 
    MoFEMFunctionReturn(0);
  };
 
  auto set_schur_pc = [&](auto solver,
                          boost::shared_ptr<SetUpSchur> &schur_ptr) {
    MoFEMFunctionBeginHot;
 
    auto name_prb = simple->getProblemName();
 
    // create sub dm for Schur complement
    auto create_schur_dm = [&](SmartPetscObj<DM> base_dm,
                               SmartPetscObj<DM> &dm_sub) {
      MoFEMFunctionBegin;
      dm_sub = createDM(mField.get_comm(), "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "SCHUR");
      CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
      CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
      CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
      for (auto f : {"U"}) {
        CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
        CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
      }
      CHKERR DMSetUp(dm_sub);
 
      MoFEMFunctionReturn(0);
    };
 
    auto create_block_dm = [&](SmartPetscObj<DM> base_dm,
                               SmartPetscObj<DM> &dm_sub) {
      MoFEMFunctionBegin;
      dm_sub = createDM(mField.get_comm(), "DMMOFEM");
      CHKERR DMMoFEMCreateSubDM(dm_sub, base_dm, "BLOCK");
      CHKERR DMMoFEMSetSquareProblem(dm_sub, PETSC_TRUE);
      CHKERR DMMoFEMAddElement(dm_sub, simple->getDomainFEName());
      CHKERR DMMoFEMAddElement(dm_sub, simple->getBoundaryFEName());
#ifdef ADD_CONTACT
      for (auto f : {"SIGMA", "EP", "TAU"}) {
        CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
        CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
      }
#else
      for (auto f : {"EP", "TAU"}) {
        CHKERR DMMoFEMAddSubFieldRow(dm_sub, f);
        CHKERR DMMoFEMAddSubFieldCol(dm_sub, f);
      }
#endif
      CHKERR DMSetUp(dm_sub);
      MoFEMFunctionReturn(0);
    };
 
    // Create nested (sub BC) Schur DM
    if constexpr (AT == AssemblyType::BLOCK_SCHUR) {
 
      SmartPetscObj<DM> dm_schur;
      CHKERR create_schur_dm(simple->getDM(), dm_schur);
      SmartPetscObj<DM> dm_block;
      CHKERR create_block_dm(simple->getDM(), dm_block);
 
#ifdef ADD_CONTACT
 
      auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
        auto block_mat_data = createBlockMatStructure(
            simple->getDM(),
 
            {
 
                {simple->getDomainFEName(),
 
                 {{"U", "U"},
                  {"SIGMA", "SIGMA"},
                  {"U", "SIGMA"},
                  {"SIGMA", "U"},
                  {"EP", "EP"},
                  {"TAU", "TAU"},
                  {"U", "EP"},
                  {"EP", "U"},
                  {"EP", "TAU"},
                  {"TAU", "EP"},
                  {"TAU", "U"}
 
                 }},
 
                {simple->getBoundaryFEName(),
 
                 {{"SIGMA", "SIGMA"}, {"U", "SIGMA"}, {"SIGMA", "U"}
 
                 }}
 
            }
 
        );
 
        return createSchurNestedMatrixStruture(
 
            {dm_schur, dm_block}, block_mat_data,
 
            {"SIGMA", "EP", "TAU"}, {nullptr, nullptr, nullptr}, true
 
        );
      };
 
#else
 
      auto get_nested_mat_data = [&](auto schur_dm, auto block_dm) {
        auto block_mat_data =
            createBlockMatStructure(simple->getDM(),
 
                                    {{simple->getDomainFEName(),
 
                                      {{"U", "U"},
                                       {"EP", "EP"},
                                       {"TAU", "TAU"},
                                       {"U", "EP"},
                                       {"EP", "U"},
                                       {"EP", "TAU"},
                                       {"TAU", "U"},
                                       {"TAU", "EP"}
 
                                      }}}
 
            );
 
        return createSchurNestedMatrixStruture(
 
            {dm_schur, dm_block}, block_mat_data,
 
            {"EP", "TAU"}, {nullptr, nullptr}, false
 
        );
      };
 
#endif
 
      auto nested_mat_data = get_nested_mat_data(dm_schur, dm_block);
      CHKERR DMMoFEMSetNestSchurData(simple->getDM(), nested_mat_data);
 
      auto block_is = getDMSubData(dm_block)->getSmartRowIs();
      auto ao_schur = getDMSubData(dm_schur)->getSmartRowMap();
 
      // Indices has to be map fro very to level, while assembling Schur
      // complement.
      schur_ptr =
          SetUpSchur::createSetUpSchur(mField, dm_schur, block_is, ao_schur);
      CHKERR schur_ptr->setUp(solver);
    }
 
    MoFEMFunctionReturnHot(0);
  };
 
  auto dm = simple->getDM();
  auto D = createDMVector(dm);
  auto DD = vectorDuplicate(D);
  uXScatter = scatter_create(D, 0);
  uYScatter = scatter_create(D, 1);
  if constexpr (SPACE_DIM == 3)
    uZScatter = scatter_create(D, 2);
 
  auto create_solver = [pip_mng]() {
    if (is_quasi_static == PETSC_TRUE)
      return pip_mng->createTSIM();
    else
      return pip_mng->createTSIM2();
  };
 
  auto solver = create_solver();
 
  auto active_pre_lhs = []() {
    MoFEMFunctionBegin;
    std::fill(PlasticOps::CommonData::activityData.begin(),
              PlasticOps::CommonData::activityData.end(), 0);
    MoFEMFunctionReturn(0);
  };
 
  auto active_post_lhs = [&]() {
    MoFEMFunctionBegin;
    auto get_iter = [&]() {
      SNES snes;
      CHK_THROW_MESSAGE(TSGetSNES(solver, &snes), "Can not get SNES");
      int iter;
      CHK_THROW_MESSAGE(SNESGetIterationNumber(snes, &iter),
                        "Can not get iter");
      return iter;
    };
 
    auto iter = get_iter();
    if (iter >= 0) {
 
      std::array<int, 5> activity_data;
      std::fill(activity_data.begin(), activity_data.end(), 0);
      MPI_Allreduce(PlasticOps::CommonData::activityData.data(),
                    activity_data.data(), activity_data.size(), MPI_INT,
                    MPI_SUM, mField.get_comm());
 
      int &active_points = activity_data[0];
      int &avtive_full_elems = activity_data[1];
      int &avtive_elems = activity_data[2];
      int &nb_points = activity_data[3];
      int &nb_elements = activity_data[4];
 
      if (nb_points) {
 
        double proc_nb_points =
            100 * static_cast<double>(active_points) / nb_points;
        double proc_nb_active =
            100 * static_cast<double>(avtive_elems) / nb_elements;
        double proc_nb_full_active = 100;
        if (avtive_elems)
          proc_nb_full_active =
              100 * static_cast<double>(avtive_full_elems) / avtive_elems;
 
        MOFEM_LOG_C("PLASTICITY", Sev::inform,
                    "Iter %d nb pts %d nb active pts %d (%3.3f\%) nb active "
                    "elements %d "
                    "(%3.3f\%) nb full active elems %d (%3.3f\%)",
                    iter, nb_points, active_points, proc_nb_points,
                    avtive_elems, proc_nb_active, avtive_full_elems,
                    proc_nb_full_active, iter);
      }
    }
 
    MoFEMFunctionReturn(0);
  };
 
  auto add_active_dofs_elem = [&](auto dm) {
    MoFEMFunctionBegin;
    auto fe_pre_proc = boost::make_shared<FEMethod>();
    fe_pre_proc->preProcessHook = active_pre_lhs;
    auto fe_post_proc = boost::make_shared<FEMethod>();
    fe_post_proc->postProcessHook = active_post_lhs;
    auto ts_ctx_ptr = getDMTsCtx(dm);
    ts_ctx_ptr->getPreProcessIJacobian().push_front(fe_pre_proc);
    ts_ctx_ptr->getPostProcessIJacobian().push_back(fe_post_proc);
    MoFEMFunctionReturn(0);
  };
 
  auto set_essential_bc = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    // This is low level pushing finite elements (pipelines) to solver
 
    auto pre_proc_ptr = boost::make_shared<FEMethod>();
    auto post_proc_rhs_ptr = boost::make_shared<FEMethod>();
    auto post_proc_lhs_ptr = boost::make_shared<FEMethod>();
 
    // Add boundary condition scaling
    auto disp_time_scale = boost::make_shared<TimeScale>();
 
    auto get_bc_hook_rhs = [this, pre_proc_ptr, disp_time_scale]() {
      EssentialPreProc<DisplacementCubitBcData> hook(mField, pre_proc_ptr,
                                                     {disp_time_scale}, false);
      return hook;
    };
    pre_proc_ptr->preProcessHook = get_bc_hook_rhs();
 
    auto get_post_proc_hook_rhs = [this, post_proc_rhs_ptr]() {
      MoFEMFunctionBegin;
      CHKERR EssentialPreProcReaction<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, nullptr, Sev::verbose)();
      CHKERR EssentialPostProcRhs<DisplacementCubitBcData>(
          mField, post_proc_rhs_ptr, 1.)();
      MoFEMFunctionReturn(0);
    };
    auto get_post_proc_hook_lhs = [this, post_proc_lhs_ptr]() {
      return EssentialPostProcLhs<DisplacementCubitBcData>(
          mField, post_proc_lhs_ptr, 1.);
    };
    post_proc_rhs_ptr->postProcessHook = get_post_proc_hook_rhs;
 
    auto ts_ctx_ptr = getDMTsCtx(dm);
    ts_ctx_ptr->getPreProcessIFunction().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPreProcessIJacobian().push_front(pre_proc_ptr);
    ts_ctx_ptr->getPostProcessIFunction().push_back(post_proc_rhs_ptr);
    post_proc_lhs_ptr->postProcessHook = get_post_proc_hook_lhs();
    ts_ctx_ptr->getPostProcessIJacobian().push_back(post_proc_lhs_ptr);
 
    MoFEMFunctionReturn(0);
  };
 
  auto B = createDMMatrix(dm);
  if (is_quasi_static == PETSC_FALSE) {
    CHKERR TSSetIJacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
  } else {
    CHKERR TSSetI2Jacobian(solver, B, B, PETSC_NULLPTR, PETSC_NULLPTR);
  }
  if (is_quasi_static == PETSC_TRUE) {
    CHKERR TSSetSolution(solver, D);
  } else {
    CHKERR TS2SetSolution(solver, D, DD);
  }
  CHKERR set_section_monitor(solver);
  CHKERR set_time_monitor(dm, solver);
  CHKERR TSSetFromOptions(solver);
 
  CHKERR add_active_dofs_elem(dm);
  boost::shared_ptr<SetUpSchur> schur_ptr;
  CHKERR set_schur_pc(solver, schur_ptr);
  CHKERR set_essential_bc(dm, solver);
 
  MOFEM_LOG_CHANNEL("TIMER");
  MOFEM_LOG_TAG("TIMER", "timer");
  if (set_timer)
    BOOST_LOG_SCOPED_THREAD_ATTR("Timeline", attrs::timer());
  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp";
  CHKERR TSSetUp(solver);
  MOFEM_LOG("TIMER", Sev::verbose) << "TSSetUp <= done";
  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve";
  CHKERR TSSolve(solver, NULL);
  MOFEM_LOG("TIMER", Sev::verbose) << "TSSolve <= done";
 
  MoFEMFunctionReturn(0);
}

Friends And Related Symbol Documentation

TSPrePostProc

friend struct TSPrePostProc

friend

Definition at line 434 of file dynamic_first_order_con_law.cpp.

Member Data Documentation

adjointDM

SmartPetscObj<DM> Example::adjointDM

private

Data manager for adjoint problem.

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 212 of file adjoint.cpp.

approxFunction

ApproxFieldFunction< FIELD_DIM > Example::approxFunction

staticprivate

Initial value:

=
    ApproxFieldFunction<FIELD_DIM>()

Examples

mofem/tutorials/scl-0/approximaton.cpp.

Definition at line 61 of file approximaton.cpp.

approxGradVals

boost::shared_ptr<MatrixDouble> Example::approxGradVals

private

Examples

mofem/tutorials/scl-8/radiation.cpp.

Definition at line 63 of file radiation.cpp.

approxVals

boost::shared_ptr<VectorDouble> Example::approxVals

private

Examples

mofem/tutorials/scl-8/radiation.cpp.

Definition at line 62 of file radiation.cpp.

aveMaxMin

std::array< double, Example::BoundingBox::LAST_BB > Example::aveMaxMin

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 110 of file phase.cpp.

base

FieldApproximationBase Example::base

private

Choice of finite element basis functions

Examples

mofem/tutorials/adv-0/plastic.cpp, mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/fun-2/plot_base.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 68 of file plot_base.cpp.

boundaryMarker

boost::shared_ptr< std::vector< unsigned char > > Example::boundaryMarker

private

Examples

mofem/tutorials/clx-0/helmholtz.cpp, and mofem/tutorials/mix-1/phase.cpp.

Definition at line 42 of file helmholtz.cpp.

commonDataPtr

boost::shared_ptr< CommonData > Example::commonDataPtr

private

Examples

mofem/tutorials/fun-1/integration.cpp, and mofem/tutorials/scl-0/approximaton.cpp.

Definition at line 42 of file integration.cpp.

domianLhsFEPtr

boost::shared_ptr<FEMethod> Example::domianLhsFEPtr

private

Examples

mofem/tutorials/vec-4/shallow_wave.cpp.

Definition at line 355 of file shallow_wave.cpp.

domianRhsFEPtr

boost::shared_ptr<FEMethod> Example::domianRhsFEPtr

private

Examples

mofem/tutorials/vec-4/shallow_wave.cpp.

Definition at line 356 of file shallow_wave.cpp.

ePS

SmartPetscObj<EPS> Example::ePS

private

Examples

mofem/tutorials/vec-1/eigen_elastic.cpp.

Definition at line 69 of file eigen_elastic.cpp.

fieldOrder

int Example::fieldOrder = 2

private

Polynomial order for approximation.

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 208 of file adjoint.cpp.

focalIndex

int Example::focalIndex

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 112 of file phase.cpp.

iI

std::vector< MatrixInt > Example::iI

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 109 of file phase.cpp.

initialGeometry

SmartPetscObj<Vec> Example::initialGeometry

private

Initial geometry field.

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 219 of file adjoint.cpp.

K

SmartPetscObj<Mat> Example::K

private

Examples

mofem/tutorials/vec-1/eigen_elastic.cpp.

Definition at line 68 of file eigen_elastic.cpp.

kspElastic

SmartPetscObj<KSP> Example::kspElastic

private

Linear solver for elastic problem.

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 211 of file adjoint.cpp.

M

SmartPetscObj<Mat> Example::M

private

Examples

mofem/tutorials/vec-1/eigen_elastic.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 67 of file eigen_elastic.cpp.

matDPtr

boost::shared_ptr<MatrixDouble> Example::matDPtr

private

Examples

mofem/tutorials/vec-1/eigen_elastic.cpp.

Definition at line 65 of file eigen_elastic.cpp.

meshVolumeAndCount

std::array<double, 2> Example::meshVolumeAndCount = {0, 0}

inlinestatic

Examples

mofem/tutorials/adv-0/plastic.cpp.

Definition at line 223 of file plastic.cpp.

{0, 0};

mField

MoFEM::Interface & Example::mField

private

Reference to MoFEM interface.

Examples

mofem/tutorials/adv-0/plastic.cpp, mofem/tutorials/adv-4/dynamic_first_order_con_law.cpp, mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-1/integration.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/mix-1/phase.cpp, mofem/tutorials/scl-0/approximaton.cpp, mofem/tutorials/scl-8/radiation.cpp, mofem/tutorials/scl-9/heat_method.cpp, mofem/tutorials/vec-1/eigen_elastic.cpp, mofem/tutorials/vec-3/nonlinear_dynamic_elastic.cpp, mofem/tutorials/vec-4/shallow_wave.cpp, mofem/tutorials/vec-5/free_surface.cpp, and mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 226 of file plastic.cpp.

modeBBoxes

std::vector<std::array<double, 6> > Example::modeBBoxes

private

Bounding boxes of optimization blocks.

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 218 of file adjoint.cpp.

modeCentroids

std::vector<std::array<double, 3> > Example::modeCentroids

private

Centroids of optimization blocks

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 217 of file adjoint.cpp.

modeVecs

std::vector<SmartPetscObj<Vec> > Example::modeVecs

private

Topology mode vectors (design variables)

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 216 of file adjoint.cpp.

pinchNodes

Range Example::pinchNodes

private

Examples

mofem/tutorials/scl-9/heat_method.cpp.

Definition at line 55 of file heat_method.cpp.

pythonPtr

boost::shared_ptr<ObjectiveFunctionData> Example::pythonPtr

private

Interface to Python objective function.

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 213 of file adjoint.cpp.

reactionFe

boost::shared_ptr<DomainEle> Example::reactionFe

private

Examples

mofem/tutorials/adv-0/plastic.cpp.

Definition at line 235 of file plastic.cpp.

rigidBodyMotion

std::array<SmartPetscObj<Vec>, 6> Example::rigidBodyMotion

private

Examples

mofem/tutorials/vec-1/eigen_elastic.cpp.

Definition at line 71 of file eigen_elastic.cpp.

rZ

std::vector< double > Example::rZ

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 108 of file phase.cpp.

savitzkyGolayNormalisation

int Example::savitzkyGolayNormalisation

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 113 of file phase.cpp.

savitzkyGolayWeights

const int * Example::savitzkyGolayWeights

staticprivate

Examples

mofem/tutorials/mix-1/phase.cpp.

Definition at line 114 of file phase.cpp.

simpleInterface

Simple * Example::simpleInterface

private

Examples

mofem/tutorials/clx-0/helmholtz.cpp, mofem/tutorials/fun-2/plot_base.cpp, mofem/tutorials/mix-1/phase.cpp, mofem/tutorials/scl-0/approximaton.cpp, and mofem/tutorials/scl-9/heat_method.cpp.

Definition at line 41 of file helmholtz.cpp.

space

FieldSpace Example::space

private

Examples

mofem/tutorials/fun-2/plot_base.cpp.

Definition at line 69 of file plot_base.cpp.

uXScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Example::uXScatter

private

Examples

mofem/tutorials/adv-0/plastic.cpp.

Definition at line 237 of file plastic.cpp.

uYScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Example::uYScatter

private

Examples

mofem/tutorials/adv-0/plastic.cpp.

Definition at line 238 of file plastic.cpp.

uZScatter

std::tuple<SmartPetscObj<Vec>, SmartPetscObj<VecScatter> > Example::uZScatter

private

Examples

mofem/tutorials/adv-0/plastic.cpp.

Definition at line 239 of file plastic.cpp.

vectorFieldPtr

boost::shared_ptr<MatrixDouble> Example::vectorFieldPtr = nullptr

private

Field values at evaluation points.

Examples

mofem/tutorials/vec-7/adjoint.cpp.

Definition at line 187 of file adjoint.cpp.

---

The documentation for this struct was generated from the following files:

• tutorials/adv-0/plastic.cpp
• tutorials/adv-4/dynamic_first_order_con_law.cpp
• tutorials/clx-0/helmholtz.cpp
• tutorials/fun-1/integration.cpp
• tutorials/fun-2/plot_base.cpp
• tutorials/mix-1/phase.cpp
• tutorials/scl-0/approximaton.cpp
• tutorials/scl-8/radiation.cpp
• tutorials/scl-9/heat_method.cpp
• tutorials/vec-1/eigen_elastic.cpp
• tutorials/vec-3/nonlinear_dynamic_elastic.cpp
• tutorials/vec-4/shallow_wave.cpp
• tutorials/vec-7/adjoint.cpp