#include <stdlib.h>
static char help[] =
"...\n\n";
public:
boost::shared_ptr<MatrixDouble> field_grad_mat)
const double area = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
const double an = 1. / std::sqrt(1 + t_field_grad(
i) * t_field_grad(
i));
locLhs(rr, cc) += (t_row_diff_base(
i) * t_col_diff_base(
i)) * an *
a -
(t_field_grad(
i) * t_col_diff_base(
i)) *
t_field_grad(
i) * an * an * an *
a;
++t_col_diff_base;
}
++t_row_diff_base;
}
++t_w;
++t_field_grad;
}
}
private:
};
public:
boost::shared_ptr<MatrixDouble> field_grad_mat)
const double area = getMeasure();
const int nb_integration_points = getGaussPts().size2();
auto t_w = getFTensor0IntegrationWeight();
for (int gg = 0; gg != nb_integration_points; gg++) {
const double a = t_w * area;
const double an = 1. / std::sqrt(1 + t_field_grad(
i) * t_field_grad(
i));
nf[rr] += (t_diff_base(
i) * t_field_grad(
i)) * an *
a;
++t_base;
++t_diff_base;
}
++t_w;
++t_field_grad;
}
}
private:
};
public:
private:
const double z) {
return sin(2 * M_PI * (x + y));
}
};
: mField(m_field) {}
}
}
}
};
}
auto get_ents_on_mesh_skin = [&]() {
CHKERR skin.find_skin(0, body_ents,
false, skin_ents);
ParallelComm *pcomm =
pcomm->filter_pstatus(skin_ents, PSTATUS_SHARED | PSTATUS_MULTISHARED,
PSTATUS_NOT, -1, &boundary_ents);
boundary_ents.merge(skin_verts);
return boundary_ents;
};
auto mark_boundary_dofs = [&](
Range &&skin_edges) {
auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
ProblemsManager::OR, skin_edges, *marker_ptr);
return marker_ptr;
};
}
auto add_domain_base_ops = [&](auto &pipeline) {
auto det_ptr = boost::make_shared<VectorDouble>();
auto jac_ptr = boost::make_shared<MatrixDouble>();
auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
};
auto add_domain_lhs_ops = [&](auto &pipeline) {
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
pipeline.push_back(
};
auto add_domain_rhs_ops = [&](auto &pipeline) {
auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
pipeline.push_back(
};
auto add_boundary_base_ops = [&](auto &pipeline) {};
auto add_lhs_base_ops = [&](auto &pipeline) {
"U", "U", [](const double, const double, const double) { return 1; }));
};
auto add_rhs_base_ops = [&](auto &pipeline) {
auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
"U", u_at_gauss_pts,
[](const double, const double, const double) { return 1; }));
};
add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
}
auto set_fieldsplit_preconditioner = [&](auto snes) {
KSP ksp;
CHKERR SNESGetKSP(snes, &ksp);
PC pc;
PetscBool is_pcfs = PETSC_FALSE;
PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
if (is_pcfs == PETSC_TRUE) {
auto name_prb =
simple->getProblemName();
int is_all_bc_size;
CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
<< "Field split block size " << is_all_bc_size;
CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
is_all_bc);
}
};
CHKERR SNESSetFromOptions(solver);
CHKERR set_fieldsplit_preconditioner(solver);
CHKERR SNESSolve(solver, global_rhs, global_solution);
SCATTER_REVERSE);
}
auto post_proc = boost::make_shared<PostProcEle>(
mField);
auto u_ptr = boost::make_shared<VectorDouble>();
post_proc->getOpPtrVector().push_back(
post_proc->getOpPtrVector().push_back(
new OpPPMap(post_proc->getPostProcMesh(),
post_proc->getMapGaussPts(),
{{"U", u_ptr}},
{}, {}, {}
)
);
CHKERR post_proc->writeFile(
"out_result.h5m");
}
int main(
int argc,
char *argv[]) {
auto core_log = logging::core::get();
core_log->add_sink(
LogManager::createSink(LogManager::getStrmWorld(), "EXAMPLE"));
LogManager::setLog("EXAMPLE");
try {
DMType dm_name = "DMMOFEM";
moab::Core mb_instance;
moab::Interface &moab = mb_instance;
CHKERR minimal_surface_problem.runProgram();
}
return 0;
}
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.
#define MOFEM_LOG(channel, severity)
Log.
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
FTensor::Index< 'i', SPACE_DIM > i
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 > smartVectorDuplicate(SmartPetscObj< 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.
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.
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
MoFEMErrorCode iNtegrate(EntData &data)
Class dedicated to integrate operator.
Integrate the domain tangent matrix (LHS)
boost::shared_ptr< MatrixDouble > fieldGradMat
MoFEMErrorCode iNtegrate(EntData &row_data, EntData &col_data)
Integrate grad-grad operator.