5static char help[] =
"...\n\n";
43 boost::shared_ptr<MatrixDouble> field_grad_mat)
54 const double area = getMeasure();
57 const int nb_integration_points = getGaussPts().size2();
59 auto t_w = getFTensor0IntegrationWeight();
61 auto t_field_grad = getFTensor1FromMat<2>(*
fieldGradMat);
67 for (
int gg = 0; gg != nb_integration_points; gg++) {
69 const double a = t_w * area;
70 const double an = 1. / std::sqrt(1 + t_field_grad(
i) * t_field_grad(
i));
79 locLhs(rr, cc) += (t_row_diff_base(
i) * t_col_diff_base(
i)) * an *
a -
81 (t_field_grad(
i) * t_col_diff_base(
i)) *
82 t_field_grad(
i) * an * an * an *
a;
121 boost::shared_ptr<MatrixDouble> field_grad_mat)
131 const double area = getMeasure();
134 const int nb_integration_points = getGaussPts().size2();
136 auto t_w = getFTensor0IntegrationWeight();
138 auto t_field_grad = getFTensor1FromMat<2>(*
fieldGradMat);
146 for (
int gg = 0; gg != nb_integration_points; gg++) {
148 const double a = t_w * area;
149 const double an = 1. / std::sqrt(1 + t_field_grad(
i) * t_field_grad(
i));
155 nf[rr] += (t_diff_base(
i) * t_field_grad(
i)) * an *
a;
196 return sin(2 * M_PI * (x + y));
271 auto get_ents_on_mesh_skin = [&]() {
276 CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
279 ParallelComm *pcomm =
281 pcomm->filter_pstatus(skin_ents, PSTATUS_SHARED | PSTATUS_MULTISHARED,
282 PSTATUS_NOT, -1, &boundary_ents);
285 mField.
get_moab().get_connectivity(boundary_ents, skin_verts,
true);
286 boundary_ents.merge(skin_verts);
288 return boundary_ents;
291 auto mark_boundary_dofs = [&](
Range &&skin_edges) {
293 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
309 auto add_domain_base_ops = [&](
auto &pipeline) {
310 auto det_ptr = boost::make_shared<VectorDouble>();
311 auto jac_ptr = boost::make_shared<MatrixDouble>();
312 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
319 auto add_domain_lhs_ops = [&](
auto &pipeline) {
321 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
329 auto add_domain_rhs_ops = [&](
auto &pipeline) {
331 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
338 auto add_boundary_base_ops = [&](
auto &pipeline) {};
340 auto add_lhs_base_ops = [&](
auto &pipeline) {
343 "U",
"U", [](
const double,
const double,
const double) {
return 1; }));
346 auto add_rhs_base_ops = [&](
auto &pipeline) {
348 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
352 [](
const double,
const double,
const double) {
return 1; }));
358 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
359 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
360 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
361 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
363 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
364 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
365 add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
366 add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
376 auto set_fieldsplit_preconditioner = [&](
auto snes) {
379 CHKERR SNESGetKSP(snes, &ksp);
381 CHKERR KSPGetPC(ksp, &pc);
382 PetscBool is_pcfs = PETSC_FALSE;
383 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
385 if (is_pcfs == PETSC_TRUE) {
387 auto name_prb =
simple->getProblemName();
392 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
394 <<
"Field split block size " << is_all_bc_size;
395 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
402 auto dm =
simple->getDM();
410 CHKERR SNESSetFromOptions(solver);
411 CHKERR set_fieldsplit_preconditioner(solver);
415 CHKERR SNESSolve(solver, global_rhs, global_solution);
427 auto post_proc = boost::make_shared<PostProcEle>(
mField);
429 auto u_ptr = boost::make_shared<VectorDouble>();
430 post_proc->getOpPtrVector().push_back(
435 post_proc->getOpPtrVector().push_back(
437 new OpPPMap(post_proc->getPostProcMesh(),
438 post_proc->getMapGaussPts(),
449 auto dm =
simple->getDM();
452 CHKERR post_proc->writeFile(
"out_result.h5m");
457int main(
int argc,
char *argv[]) {
464 auto core_log = logging::core::get();
473 DMType dm_name =
"DMMOFEM";
477 moab::Core mb_instance;
478 moab::Interface &moab = mb_instance;
void simple(double P1[], double P2[], double P3[], double c[], const int N)
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
ElementsAndOps< SPACE_DIM >::BoundaryEle BoundaryEle
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MYPCOMM_INDEX
default communicator number PCOMM
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
PetscErrorCode DMoFEMMeshToGlobalVector(DM dm, Vec g, InsertMode mode, ScatterMode scatter_mode)
set ghosted vector values on all existing mesh entities
SmartPetscObj< SNES > createSNES(SmartPetscObj< DM > dm=nullptr)
Create SNES (nonlinear) solver.
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FTensor::Index< 'i', 2 > i
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
implementation of Data Operators for Forces and Sources
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< G >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
constexpr auto field_name
static constexpr int approx_order
OpBaseImpl< PETSC, EdgeEleOp > OpBase
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode outputResults()
MoFEMErrorCode setIntegrationRules()
MoFEMErrorCode runProgram()
MoFEMErrorCode solveSystem()
MoFEM::Interface & mField
static double boundaryFunction(const double x, const double y, const double z)
MinimalSurfaceEqn(MoFEM::Interface &m_field)
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
MoFEMErrorCode setupProblem()
Simple interface for fast problem set-up.
virtual moab::Interface & get_moab()=0
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
FTensor::Tensor1< FTensor::PackPtr< double *, Tensor_Dim >, Tensor_Dim > getFTensor1DiffN(const FieldApproximationBase base)
Get derivatives of base functions.
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
default operator for TRI element
Section manager is used to create indexes and sections.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
static boost::shared_ptr< std::ostream > getStrmWorld()
Get the strm world object.
VectorDouble locF
local entity vector
int nbRows
number of dofs on rows
MatrixDouble locMat
local entity block matrix
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
MoFEM::FaceElementForcesAndSourcesCore FaceEle
MoFEM::EdgeElementForcesAndSourcesCore EdgeEle
Problem manager is used to build and partition problems.
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Simple interface for fast problem set-up.
intrusive_ptr for managing petsc objects
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Integrate the domain residual vector (RHS)
boost::shared_ptr< MatrixDouble > fieldGradMat
OpDomainResidualVector(std::string field_name, boost::shared_ptr< MatrixDouble > field_grad_mat)
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
Integrate the domain tangent matrix (LHS)
boost::shared_ptr< MatrixDouble > fieldGradMat
OpDomainTangentMatrix(std::string row_field_name, std::string col_field_name, boost::shared_ptr< MatrixDouble > field_grad_mat)
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.