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 Member Functions
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	[Run problem] More...
	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 ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()
	Example (MoFEM::Interface &m_field)
MoFEMErrorCode	runProblem ()

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] More...
MoFEMErrorCode	createCommonData ()
	[Set up problem] More...
MoFEMErrorCode	bC ()
	[Create common data] More...
MoFEMErrorCode	OPs ()
	[Boundary condition] More...
MoFEMErrorCode	tsSolve ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	bC ()
MoFEMErrorCode	OPs ()
MoFEMErrorCode	tsSolve ()
MoFEMErrorCode	readMesh ()
	[Run problem] More...
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	boundaryCondition ()
	[Set up problem] More...
MoFEMErrorCode	assembleSystem ()
	[Push operators to pipeline] More...
MoFEMErrorCode	solveSystem ()
	[Solve] More...
MoFEMErrorCode	outputResults ()
	[Solve] More...
MoFEMErrorCode	checkResults ()
	[Postprocess results] More...
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	setUp ()
	[Run all] More...
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	setFieldValues ()
	[Create common data] More...
MoFEMErrorCode	pushOperators ()
	[Set density distribution] More...
MoFEMErrorCode	integrateElements ()
	[Push operators to pipeline] More...
MoFEMErrorCode	postProcess ()
	[Integrate] More...
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	setIntegrationRules ()
	[Set up problem] More...
MoFEMErrorCode	createCommonData ()
MoFEMErrorCode	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	readMesh ()
MoFEMErrorCode	setupProblem ()
MoFEMErrorCode	assembleSystemIntensity ()
	[Calculate flux on boundary] More...
MoFEMErrorCode	assembleSystemFlux ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	calculateFlux (double &calc_flux)
	[Set up problem] More...
MoFEMErrorCode	outputResults (const int i)
	[Solve] More...
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] More...
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	boundaryCondition ()
MoFEMErrorCode	assembleSystem ()
MoFEMErrorCode	solveSystem ()
MoFEMErrorCode	outputResults ()
MoFEMErrorCode	checkResults ()
MoFEMErrorCode	addMatBlockOps (boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &pipeline, std::string field_name, std::string block_name, boost::shared_ptr< MatrixDouble > mat_D_Ptr, Sev sev)
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	gettingNorms ()
	[Solve] More...
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 ()

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
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
FieldSpace	space
boost::shared_ptr< VectorDouble >	approxVals
boost::shared_ptr< MatrixDouble >	approxGradVals
Range	pinchNodes
PetscBool	doEvalField
std::array< double, SPACE_DIM >	fieldEvalCoords
boost::shared_ptr< FieldEvaluatorInterface::SetPtsData >	fieldEvalData
boost::shared_ptr< MatrixDouble >	vectorFieldPtr
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

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]

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 177 of file plastic.cpp.

Member Enumeration Documentation

BoundingBox

enum Example::BoundingBox

private

Enumerator
CENTER_X
CENTER_Y
MAX_X
MAX_Y
MIN_X
MIN_Y
LAST_BB

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/15]

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

inline

Definition at line 179 of file plastic.cpp.

: mField(m_field) {}

Example() [2/15]

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

inline

Definition at line 230 of file plasticity_old_schur.cpp.

: mField(m_field) {}

Example() [3/15]

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

inline

Definition at line 421 of file dynamic_first_order_con_law.cpp.

: mField(m_field) {}

Example() [4/15]

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

inline

Definition at line 35 of file helmholtz.cpp.

: mField(m_field) {}

Example() [5/15]

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

inline

Definition at line 25 of file integration.cpp.

: mField(m_field) {}

Example() [6/15]

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

inline

Definition at line 50 of file plot_base.cpp.

: mField(m_field) {}

Example() [7/15]

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

inline

Definition at line 81 of file phase.cpp.

: mField(m_field) {}

Example() [8/15]

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

inline

Definition at line 53 of file approximaton.cpp.

: mField(m_field) {}

Example() [9/15]

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

inline

Definition at line 45 of file radiation.cpp.

: mField(m_field) {}

Example() [10/15]

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

inline

Definition at line 46 of file heat_method.cpp.

: mField(m_field) {}

Example() [11/15]

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

inline

Definition at line 88 of file elastic.cpp.

: mField(m_field) {}

Example() [12/15]

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

inline

Definition at line 49 of file eigen_elastic.cpp.

: mField(m_field) {}

Example() [13/15]

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

inline

Definition at line 59 of file nonlinear_elastic.cpp.

: mField(m_field) {}

Example() [14/15]

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

inline

Definition at line 55 of file nonlinear_dynamic_elastic.cpp.

: mField(m_field) {}

Example() [15/15]

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

inline

Definition at line 339 of file shallow_wave.cpp.

: mField(m_field) {}

Member Function Documentation

addMatBlockOps()

MoFEMErrorCode Example::addMatBlockOps ( boost::ptr_deque< ForcesAndSourcesCore::UserDataOperator > &  pipeline, std::string  field_name, std::string  block_name, boost::shared_ptr< MatrixDouble >  mat_D_Ptr, Sev  sev  )

private

[Calculate elasticity tensor]

[Calculate elasticity tensor]

Definition at line 114 of file elastic.cpp.

                                                      {
  MoFEMFunctionBegin;
 
  struct OpMatBlocks : public DomainEleOp {
    OpMatBlocks(std::string field_name, boost::shared_ptr<MatrixDouble> m,
                double bulk_modulus_K, double shear_modulus_G,
                MoFEM::Interface &m_field, Sev sev,
                std::vector<const CubitMeshSets *> meshset_vec_ptr)
        : DomainEleOp(field_name, DomainEleOp::OPROW), matDPtr(m),
          bulkModulusKDefault(bulk_modulus_K),
          shearModulusGDefault(shear_modulus_G) {
      std::fill(&(doEntities[MBEDGE]), &(doEntities[MBMAXTYPE]), false);
      CHK_THROW_MESSAGE(extractBlockData(m_field, meshset_vec_ptr, sev),
                        "Can not get data from block");
    }
 
    MoFEMErrorCode doWork(int side, EntityType type,
                          EntitiesFieldData::EntData &data) {
      MoFEMFunctionBegin;
 
      for (auto &b : blockData) {
 
        if (b.blockEnts.find(getFEEntityHandle()) != b.blockEnts.end()) {
          CHKERR getMatDPtr(matDPtr, b.bulkModulusK, b.shearModulusG);
          MoFEMFunctionReturnHot(0);
        }
      }
 
      CHKERR getMatDPtr(matDPtr, bulkModulusKDefault, shearModulusGDefault);
      MoFEMFunctionReturn(0);
    }
 
  private:
    boost::shared_ptr<MatrixDouble> matDPtr;
 
    struct BlockData {
      double bulkModulusK;
      double shearModulusG;
      Range blockEnts;
    };
 
    double bulkModulusKDefault;
    double shearModulusGDefault;
    std::vector<BlockData> blockData;
 
    MoFEMErrorCode
    extractBlockData(MoFEM::Interface &m_field,
                     std::vector<const CubitMeshSets *> meshset_vec_ptr,
                     Sev sev) {
      MoFEMFunctionBegin;
 
      for (auto m : meshset_vec_ptr) {
        MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock") << *m;
        std::vector<double> block_data;
        CHKERR m->getAttributes(block_data);
        if (block_data.size() < 2) {
          SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
                  "Expected that block has two attributes");
        }
        auto get_block_ents = [&]() {
          Range ents;
          CHKERR
          m_field.get_moab().get_entities_by_handle(m->meshset, ents, true);
          return ents;
        };
 
        double young_modulus = block_data[0];
        double poisson_ratio = block_data[1];
        double bulk_modulus_K = young_modulus / (3 * (1 - 2 * poisson_ratio));
        double shear_modulus_G = young_modulus / (2 * (1 + poisson_ratio));
 
        MOFEM_TAG_AND_LOG("WORLD", sev, "MatBlock")
            << "E = " << young_modulus << " nu = " << poisson_ratio;
 
        blockData.push_back(
            {bulk_modulus_K, shear_modulus_G, get_block_ents()});
      }
      MOFEM_LOG_CHANNEL("WORLD");
      MoFEMFunctionReturn(0);
    }
 
    MoFEMErrorCode getMatDPtr(boost::shared_ptr<MatrixDouble> mat_D_ptr,
                              double bulk_modulus_K, double shear_modulus_G) {
      MoFEMFunctionBegin;
      //! [Calculate elasticity tensor]
      auto set_material_stiffness = [&]() {
        FTensor::Index<'i', SPACE_DIM> i;
        FTensor::Index<'j', SPACE_DIM> j;
        FTensor::Index<'k', SPACE_DIM> k;
        FTensor::Index<'l', SPACE_DIM> l;
        constexpr auto t_kd = FTensor::Kronecker_Delta_symmetric<int>();
        double A = 1.;
        if (SPACE_DIM == 2 && !is_plane_strain) {
          A = 2 * shear_modulus_G /
              (bulk_modulus_K + (4. / 3.) * shear_modulus_G);
        }
        auto t_D = getFTensor4DdgFromMat<SPACE_DIM, SPACE_DIM, 0>(*mat_D_ptr);
        t_D(i, j, k, l) =
            2 * shear_modulus_G * ((t_kd(i, k) ^ t_kd(j, l)) / 4.) +
            A * (bulk_modulus_K - (2. / 3.) * shear_modulus_G) * t_kd(i, j) *
                t_kd(k, l);
      };
      //! [Calculate elasticity tensor]
      constexpr auto size_symm = (SPACE_DIM * (SPACE_DIM + 1)) / 2;
      mat_D_ptr->resize(size_symm * size_symm, 1);
      set_material_stiffness();
      MoFEMFunctionReturn(0);
    }
  };
 
  pipeline.push_back(new OpMatBlocks(
      field_name, mat_D_Ptr, bulk_modulus_K, shear_modulus_G, mField, sev,
 
      // Get blockset using regular expression
      mField.getInterface<MeshsetsManager>()->getCubitMeshsetPtr(std::regex(
 
          (boost::format("%s(.*)") % block_name).str()
 
              ))
 
          ));
 
  MoFEMFunctionReturn(0);
}

assembleSystem() [1/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [2/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [3/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [4/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [5/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [6/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [7/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [8/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [9/10]

MoFEMErrorCode Example::assembleSystem ( )

private

assembleSystem() [10/10]

MoFEMErrorCode Example::assembleSystem ( )

private

[Push operators to pipeline]

[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 Internal forces]

[Push Internal forces]

[Push Body forces]

[Push Body forces]

[Push natural boundary conditions]

[Push natural boundary conditions]

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 647 of file dynamic_first_order_con_law.cpp.

                                       {
  MoFEMFunctionBegin;
  auto *simple = mField.getInterface<Simple>();
  auto *pipeline_mng = mField.getInterface<PipelineManager>();
 
  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_NULL, "", "-unit_weight", &unit_weight,
                                 PETSC_NULL);
      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_NULL, "", "-tau", &tau, PETSC_NULL);
 
    double xi = 0.;
    CHKERR PetscOptionsGetReal(PETSC_NULL, "", "-xi", &xi, PETSC_NULL);
 
    // 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_NULL, "", "-is_linear_elasticity",
                               &is_linear_elasticity, PETSC_NULL);
 
    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);
  };
 
  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);
}

assembleSystemFlux()

MoFEMErrorCode Example::assembleSystemFlux ( )

private

assembleSystemIntensity()

MoFEMErrorCode Example::assembleSystemIntensity ( )

private

[Calculate flux on boundary]

[Push operators to pipeline]

Examples

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/3]

MoFEMErrorCode Example::bC ( )

private

bC() [2/3]

MoFEMErrorCode Example::bC ( )

private

[Create common data]

[Boundary condition]

Examples

plastic.cpp.

Definition at line 504 of file plastic.cpp.

                           {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
  auto prb_mng = mField.getInterface<ProblemsManager>();
 
  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() [3/3]

MoFEMErrorCode Example::bC ( )

private

boundaryCondition() [1/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [2/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [3/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [4/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [5/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [6/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [7/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [8/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [9/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

boundaryCondition() [10/10]

MoFEMErrorCode Example::boundaryCondition ( )

private

[Set up problem]

[Create common data]

[Boundary condition]

[Applying essential BC]

[Define gravity vector]

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 536 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_NULL, "", "-sin_time_function",
                             &sin_time_function, PETSC_NULL);
 
  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);
}

calculateFlux()

MoFEMErrorCode Example::calculateFlux ( double &  calc_flux )

private

[Set up problem]

[Calculate flux on boundary]

Examples

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);
}

checkResults() [1/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [2/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [3/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [4/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [5/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [6/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [7/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [8/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [9/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [10/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [11/12]

MoFEMErrorCode Example::checkResults ( )

private

checkResults() [12/12]

MoFEMErrorCode Example::checkResults ( )

private

[Postprocess results]

[Postprocessing results]

[Print results]

[Check]

[Check results]

[Test example]

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 1205 of file dynamic_first_order_con_law.cpp.

                                     {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

createCommonData() [1/9]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [2/9]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [3/9]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [4/9]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [5/9]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [6/9]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [7/9]

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

eigen_elastic.cpp, heat_method.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 396 of file plastic.cpp.

                                         {
  MoFEMFunctionBegin;
 
  auto get_command_line_parameters = [&]() {
    MoFEMFunctionBegin;
 
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-scale", &scale, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-young_modulus",
                                 &young_modulus, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-poisson_ratio",
                                 &poisson_ratio, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening", &H, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-hardening_viscous", &visH,
                                 PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-yield_stress", &sigmaY,
                                 PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn0", &cn0, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn1", &cn1, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-zeta", &zeta, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-Qinf", &Qinf, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-b_iso", &b_iso, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-C1_k", &C1_k, PETSC_NULL);
    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-large_strains",
                               &is_large_strains, PETSC_NULL);
    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-set_timer", &set_timer,
                               PETSC_NULL);
 
    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
    PetscBool tau_order_is_set; ///< true if tau order is set
    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-tau_order", &tau_order,
                              &tau_order_is_set);
    PetscBool ep_order_is_set; ///< true if tau order is set
    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-ep_order", &ep_order,
                              &ep_order_is_set);
    CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-geom_order", &geom_order,
                              PETSC_NULL);
 
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-rho", &rho, PETSC_NULL);
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-alpha_damping",
                                 &alpha_damping, PETSC_NULL);
 
    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_NULL, "", "-is_scale", &is_scale,
                               PETSC_NULL);
    if (is_scale) {
      scale /= young_modulus;
    }
 
    MOFEM_LOG("PLASTICITY", Sev::inform) << "Scale " << scale;
 
#ifdef ADD_CONTACT
    CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-cn_contact",
                                 &ContactOps::cn_contact, PETSC_NULL);
    MOFEM_LOG("CONTACT", Sev::inform)
        << "cn_contact " << ContactOps::cn_contact;
#endif // ADD_CONTACT
 
    CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-quasi_static",
                               &is_quasi_static, PETSC_NULL);
    MOFEM_LOG("PLASTICITY", Sev::inform)
        << "Is quasi static: " << (is_quasi_static ? "true" : "false");
 
    MoFEMFunctionReturn(0);
  };
 
  CHKERR get_command_line_parameters();
 
#ifdef ADD_CONTACT
#ifdef PYTHON_SDF
  sdfPythonPtr = boost::make_shared<ContactOps::SDFPython>();
  CHKERR sdfPythonPtr->sdfInit("sdf.py");
  ContactOps::sdfPythonWeakPtr = sdfPythonPtr;
#endif
#endif // ADD_CONTACT
 
  MoFEMFunctionReturn(0);
}

createCommonData() [8/9]

MoFEMErrorCode Example::createCommonData ( )

private

createCommonData() [9/9]

MoFEMErrorCode Example::createCommonData ( )

private

getCoordsInImage()

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

staticprivate

Examples

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;
}

gettingNorms()

MoFEMErrorCode Example::gettingNorms ( )

private

[Solve]

[Getting norms]

Examples

nonlinear_elastic.cpp.

Definition at line 383 of file nonlinear_elastic.cpp.

                                     {
  MoFEMFunctionBegin;
 
  auto simple = mField.getInterface<Simple>();
  auto dm = simple->getDM();
 
  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) << "Solution norm " << nrm2;
 
  auto post_proc_norm_fe = boost::make_shared<DomainEle>(mField);
 
  auto post_proc_norm_rule_hook = [](int, int, int p) -> int { return 2 * p; };
  post_proc_norm_fe->getRuleHook = post_proc_norm_rule_hook;
 
  CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(
      post_proc_norm_fe->getOpPtrVector(), {H1});
 
  enum NORMS { U_NORM_L2 = 0, PIOLA_NORM, LAST_NORM };
  auto norms_vec =
      createVectorMPI(mField.get_comm(),
                      (mField.get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
  CHKERR VecZeroEntries(norms_vec);
 
  auto u_ptr = boost::make_shared<MatrixDouble>();
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalculateVectorFieldValues<SPACE_DIM>("U", u_ptr));
 
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor1<SPACE_DIM>(u_ptr, norms_vec, U_NORM_L2));
 
  auto common_ptr = commonDataFactory<SPACE_DIM, GAUSS, DomainEleOp>(
      mField, post_proc_norm_fe->getOpPtrVector(), "U", "MAT_ELASTIC",
      Sev::inform);
 
  post_proc_norm_fe->getOpPtrVector().push_back(
      new OpCalcNormL2Tensor2<SPACE_DIM, SPACE_DIM>(
          common_ptr->getMatFirstPiolaStress(), norms_vec, PIOLA_NORM));
 
  CHKERR DMoFEMLoopFiniteElements(dm, simple->getDomainFEName(),
                                  post_proc_norm_fe);
 
  CHKERR VecAssemblyBegin(norms_vec);
  CHKERR VecAssemblyEnd(norms_vec);
 
  MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
  if (mField.get_comm_rank() == 0) {
    const double *norms;
    CHKERR VecGetArrayRead(norms_vec, &norms);
    MOFEM_TAG_AND_LOG("SELF", Sev::inform, "example")
        << "norm_u: " << std::scientific << std::sqrt(norms[U_NORM_L2]);
    MOFEM_TAG_AND_LOG("SELF", Sev::inform, "example")
        << "norm_piola: " << std::scientific << std::sqrt(norms[PIOLA_NORM]);
    CHKERR VecRestoreArrayRead(norms_vec, &norms);
  }
 
  MoFEMFunctionReturn(0);
}

integrateElements()

MoFEMErrorCode Example::integrateElements ( )

private

[Push operators to pipeline]

[Integrate]

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

Definition at line 52 of file radiation.cpp.

{ return 2 * p_data; };

kspSolve()

MoFEMErrorCode Example::kspSolve ( )

private

[Push operators to pipeline]

[Solve]

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

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/3]

MoFEMErrorCode Example::OPs ( )

private

OPs() [2/3]

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

plastic.cpp.

Definition at line 549 of file plastic.cpp.

                            {
  MoFEMFunctionBegin;
  auto pip_mng = mField.getInterface<PipelineManager>();
  auto simple = mField.getInterface<Simple>();
  auto bc_mng = mField.getInterface<BcManager>();
 
  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 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
    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 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 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 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 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() [3/3]

MoFEMErrorCode Example::OPs ( )

private

outputResults() [1/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [2/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [3/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [4/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [5/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [6/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [7/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [8/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [9/11]

MoFEMErrorCode Example::outputResults ( )

private

outputResults() [10/11]

MoFEMErrorCode Example::outputResults ( )

private

[Solve]

[Getting norms]

[Postprocess results]

[Postprocessing results]

[Postprocess clean]

[Postprocess clean]

[Postprocess initialise]

[Postprocess initialise]

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 1183 of file dynamic_first_order_con_law.cpp.

                                      {
  MoFEMFunctionBegin;
  PetscBool test_flg = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
  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() [11/11]

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]

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

pushOperators()

MoFEMErrorCode Example::pushOperators ( )

private

[Set density distribution]

[Push operators to pipeline]

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/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [2/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [3/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [4/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [5/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [6/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [7/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [8/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [9/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [10/11]

MoFEMErrorCode Example::readMesh ( )

private

readMesh() [11/11]

MoFEMErrorCode Example::readMesh ( )

private

[Run problem]

[Run programme]

[run problem]

[Read mesh]

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 463 of file dynamic_first_order_con_law.cpp.

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

rhsSource()

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

staticprivate

Examples

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/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [2/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [3/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [4/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [5/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [6/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [7/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [8/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [9/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [10/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [11/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [12/15]

MoFEMErrorCode Example::runProblem

(

)

[Run problem]

[Create common data]

[Run programme]

[Operator]

[run problem]

[Run all]

[Run problem]

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 213 of file plastic.cpp.

                                   {
  MoFEMFunctionBegin;
  CHKERR createCommonData();
  CHKERR setupProblem();
  CHKERR bC();
  CHKERR OPs();
  CHKERR tsSolve();
  MoFEMFunctionReturn(0);
}

runProblem() [13/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [14/15]

MoFEMErrorCode Example::runProblem

(

)

runProblem() [15/15]

MoFEMErrorCode Example::runProblem

(

)

setFieldValues()

MoFEMErrorCode Example::setFieldValues ( )

private

[Create common data]

[Set density distribution]

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

setIntegrationRules() [2/3]

MoFEMErrorCode Example::setIntegrationRules ( )

private

[Set up problem]

[Set integration rule]

Examples

heat_method.cpp, and plot_base.cpp.

Definition at line 203 of file plot_base.cpp.

                                            {
  MoFEMFunctionBegin;
  MoFEMFunctionReturn(0);
}

setIntegrationRules() [3/3]

MoFEMErrorCode Example::setIntegrationRules ( )

private

setUp()

MoFEMErrorCode Example::setUp ( )

private

[Run all]

[Set up problem]

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);
}

setupProblem() [1/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [2/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [3/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [4/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [5/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [6/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [7/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [8/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [9/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [10/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [11/14]

MoFEMErrorCode Example::setupProblem ( )

private

[Run problem]

[Read mesh]

[Set up problem]

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plastic.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 225 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_NULL, "", "-order_edge", &order_edge,
                             PETSC_NULL);
  PetscBool order_face = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_face", &order_face,
                             PETSC_NULL);
  PetscBool order_volume = PETSC_FALSE;
  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-order_volume", &order_volume,
                             PETSC_NULL);
 
  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 = false;
    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_NULL, "", "-project_geometry",
                             &project_geometry, PETSC_NULL);
  if (project_geometry) {
    CHKERR project_ho_geometry();
  }
 
  MoFEMFunctionReturn(0);
}

setupProblem() [12/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [13/14]

MoFEMErrorCode Example::setupProblem ( )

private

setupProblem() [14/14]

MoFEMErrorCode Example::setupProblem ( )

private

solveSystem() [1/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [2/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [3/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [4/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [5/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [6/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [7/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [8/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [9/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [10/11]

MoFEMErrorCode Example::solveSystem ( )

private

solveSystem() [11/11]

MoFEMErrorCode Example::solveSystem ( )

private

[Solve]

[Push operators to pipeline]

[Solve]

< Mass matrix

< Linear solver

Examples

dynamic_first_order_con_law.cpp, eigen_elastic.cpp, heat_method.cpp, helmholtz.cpp, nonlinear_elastic.cpp, phase.cpp, plot_base.cpp, and shallow_wave.cpp.

Definition at line 893 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_NULL, "", "-is_linear_elasticity",
                               &is_linear_elasticity, PETSC_NULL);
 
    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_NULL, "", "-density", &rho, PETSC_NULL);
  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>();
    if (SPACE_DIM == 3) {
      CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree3D(
          field_eval_data, simple->getDomainFEName());
    } else {
      CHKERR mField.getInterface<FieldEvaluatorInterface>()->buildTree2D(
          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);
 
  auto fb = mField.getInterface<FieldBlas>();
 
  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);
}

tsSolve() [1/2]

MoFEMErrorCode Example::tsSolve ( )

private

Examples

plastic.cpp.

Definition at line 734 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 pp_fe = boost::make_shared<PostProcEle>(mField);
    auto &pip = pp_fe->getOpPtrVector();
 
    auto push_vol_ops = [this](auto &pip) {
      CHKERR AddHOOps<SPACE_DIM, SPACE_DIM, SPACE_DIM>::add(pip, {H1, HDIV},
                                                            "GEOMETRY");
 
      auto [common_plastic_ptr, common_henky_ptr] =
          PlasticOps::createCommonPlasticOps<SPACE_DIM, IT, DomainEleOp>(
              mField, "MAT_PLASTIC", pip, "U", "EP", "TAU", 1., Sev::inform);
 
      if (common_henky_ptr) {
        if (common_plastic_ptr->mGradPtr != common_henky_ptr->matGradPtr)
          CHK_THROW_MESSAGE(MOFEM_DATA_INCONSISTENCY, "Wrong pointer for grad");
      }
 
      return std::make_pair(common_plastic_ptr, common_henky_ptr);
    };
 
    auto push_vol_post_proc_ops = [this](auto &pp_fe, auto &&p) {
      MoFEMFunctionBegin;
 
      auto &pip = pp_fe->getOpPtrVector();
 
      auto [common_plastic_ptr, common_henky_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_henky_ptr->matGradPtr},
                 {"FIRST_PIOLA", common_henky_ptr->getMatFirstPiolaStress()}},
 
                {{"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);
    };
 
    auto vol_post_proc = [this, push_vol_post_proc_ops, push_vol_ops]() {
      PetscBool post_proc_vol = PETSC_FALSE;
      CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_vol",
                                 &post_proc_vol, PETSC_NULL);
      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]() {
      PetscBool post_proc_skin = PETSC_TRUE;
      CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-post_proc_skin",
                                 &post_proc_skin, PETSC_NULL);
      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_NULL, &scatter);
    return std::make_tuple(SmartPetscObj<Vec>(v),
                           SmartPetscObj<VecScatter>(scatter));
  };
 
  auto set_time_monitor = [&](auto dm, auto solver) {
    MoFEMFunctionBegin;
    boost::shared_ptr<Monitor> monitor_ptr(
        new Monitor(dm, create_post_process_elements(), reactionFe, uXScatter,
                    uYScatter, uZScatter));
    boost::shared_ptr<ForcesAndSourcesCore> null;
    CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
                               monitor_ptr, null, null);
    MoFEMFunctionReturn(0);
  };
 
  auto set_schur_pc = [&](auto solver,
                          boost::shared_ptr<SetUpSchur> &schur_ptr) {
    MoFEMFunctionBeginHot;
 
    auto bc_mng = mField.getInterface<BcManager>();
    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 getNestSchurData(
 
            {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 getNestSchurData(
 
            {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 avtive 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);
  };
 
  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 TSSetUp(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);
}

tsSolve() [2/2]

MoFEMErrorCode Example::tsSolve ( )

private

Friends And Related Function Documentation

TSPrePostProc

friend struct TSPrePostProc

friend

Definition at line 435 of file dynamic_first_order_con_law.cpp.

Member Data Documentation

approxFunction

ApproxFieldFunction< FIELD_DIM > Example::approxFunction

staticprivate

Initial value:

=
    ApproxFieldFunction<FIELD_DIM>()

Definition at line 61 of file approximaton.cpp.

approxGradVals

boost::shared_ptr<MatrixDouble> Example::approxGradVals

private

Definition at line 63 of file radiation.cpp.

approxVals

boost::shared_ptr<VectorDouble> Example::approxVals

private

Definition at line 62 of file radiation.cpp.

aveMaxMin

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

staticprivate

Examples

phase.cpp.

Definition at line 110 of file phase.cpp.

base

FieldApproximationBase Example::base

private

Definition at line 68 of file plot_base.cpp.

boundaryMarker

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

private

Definition at line 42 of file helmholtz.cpp.

commonDataPtr

boost::shared_ptr< CommonData > Example::commonDataPtr

private

Definition at line 40 of file integration.cpp.

doEvalField

PetscBool Example::doEvalField

private

Definition at line 95 of file elastic.cpp.

domianLhsFEPtr

boost::shared_ptr<FEMethod> Example::domianLhsFEPtr

private

Definition at line 355 of file shallow_wave.cpp.

domianRhsFEPtr

boost::shared_ptr<FEMethod> Example::domianRhsFEPtr

private

Definition at line 356 of file shallow_wave.cpp.

ePS

SmartPetscObj<EPS> Example::ePS

private

Examples

eigen_elastic.cpp.

Definition at line 69 of file eigen_elastic.cpp.

fieldEvalCoords

std::array<double, SPACE_DIM> Example::fieldEvalCoords

private

Definition at line 96 of file elastic.cpp.

fieldEvalData

boost::shared_ptr<FieldEvaluatorInterface::SetPtsData> Example::fieldEvalData

private

Definition at line 97 of file elastic.cpp.

focalIndex

int Example::focalIndex

staticprivate

Examples

phase.cpp.

Definition at line 112 of file phase.cpp.

iI

std::vector< MatrixInt > Example::iI

staticprivate

Examples

phase.cpp.

Definition at line 109 of file phase.cpp.

K

SmartPetscObj<Mat> Example::K

private

Definition at line 68 of file eigen_elastic.cpp.

M

SmartPetscObj<Mat> Example::M

private

Definition at line 67 of file eigen_elastic.cpp.

matDPtr

boost::shared_ptr<MatrixDouble> Example::matDPtr

private

Definition at line 65 of file eigen_elastic.cpp.

mField

MoFEM::Interface & Example::mField

private

Examples

dynamic_first_order_con_law.cpp, free_surface.cpp, and helmholtz.cpp.

Definition at line 184 of file plastic.cpp.

pinchNodes

Range Example::pinchNodes

private

Definition at line 54 of file heat_method.cpp.

reactionFe

boost::shared_ptr< DomainEle > Example::reactionFe

private

Definition at line 192 of file plastic.cpp.

rigidBodyMotion

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

private

Definition at line 71 of file eigen_elastic.cpp.

rZ

std::vector< double > Example::rZ

staticprivate

Examples

phase.cpp.

Definition at line 108 of file phase.cpp.

savitzkyGolayNormalisation

int Example::savitzkyGolayNormalisation

staticprivate

Examples

phase.cpp.

Definition at line 113 of file phase.cpp.

savitzkyGolayWeights

const int * Example::savitzkyGolayWeights

staticprivate

Examples

phase.cpp.

Definition at line 114 of file phase.cpp.

simpleInterface

Simple * Example::simpleInterface

private

Examples

helmholtz.cpp.

Definition at line 41 of file helmholtz.cpp.

space

FieldSpace Example::space

private

Definition at line 69 of file plot_base.cpp.

uXScatter

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

private

Definition at line 194 of file plastic.cpp.

uYScatter

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

private

Definition at line 195 of file plastic.cpp.

uZScatter

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

private

Definition at line 196 of file plastic.cpp.

vectorFieldPtr

boost::shared_ptr<MatrixDouble> Example::vectorFieldPtr

private

Definition at line 98 of file elastic.cpp.

---

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

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