25 using namespace MoFEM;
27 static char help[] =
"...\n\n";
29 #include <BasicFiniteElements.hpp>
54 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
71 CHKERR boundaryCondition();
84 CHKERR mField.getInterface(simpleInterface);
85 CHKERR simpleInterface->getOptions();
86 CHKERR simpleInterface->loadFile();
96 CHKERR simpleInterface->addDomainField(
"P_REAL",
H1,
98 CHKERR simpleInterface->addDomainField(
"P_IMAG",
H1,
100 CHKERR simpleInterface->addBoundaryField(
"P_REAL",
H1,
102 CHKERR simpleInterface->addBoundaryField(
"P_IMAG",
H1,
106 CHKERR simpleInterface->setFieldOrder(
"P_REAL",
order);
107 CHKERR simpleInterface->setFieldOrder(
"P_IMAG",
order);
108 CHKERR simpleInterface->setUp();
117 auto get_ents_on_mesh_skin = [&]() {
118 Range boundary_entities;
120 std::string entity_name = it->getName();
121 if (entity_name.compare(0, 2,
"BC") == 0) {
122 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
123 boundary_entities,
true);
127 Range boundary_vertices;
128 CHKERR mField.get_moab().get_connectivity(boundary_entities,
129 boundary_vertices,
true);
130 boundary_entities.merge(boundary_vertices);
132 return boundary_entities;
135 auto mark_boundary_dofs = [&](Range &&skin_edges) {
137 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
138 problem_manager->
markDofs(simpleInterface->getProblemName(),
ROW,
139 skin_edges, *marker_ptr);
143 auto remove_dofs_from_problem = [&](Range &&skin_edges) {
146 CHKERR problem_manager->removeDofsOnEntities(
147 simpleInterface->getProblemName(),
"P_IMAG", skin_edges, 0, 1);
154 boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
155 CHKERR remove_dofs_from_problem(get_ents_on_mesh_skin());
169 auto beta = [](
const double,
const double,
const double) {
return -1; };
170 auto k2 = [
k](
const double,
const double,
const double) {
return pow(
k, 2); };
171 auto kp = [
k](
const double,
const double,
const double) {
return k; };
172 auto km = [
k](
const double,
const double,
const double) {
return -
k; };
175 auto set_domain = [&]() {
177 auto det_ptr = boost::make_shared<VectorDouble>();
178 auto jac_ptr = boost::make_shared<MatrixDouble>();
179 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
188 new OpSetBc(
"P_REAL",
true, boundaryMarker));
206 auto set_boundary = [&]() {
209 new OpSetBc(
"P_REAL",
true, boundaryMarker));
217 new OpSetBc(
"P_REAL",
false, boundaryMarker));
223 new OpSetBc(
"P_REAL",
false, boundaryMarker));
245 CHKERR KSPSetFromOptions(solver);
248 auto dm = simpleInterface->getDM();
253 CHKERR VecGhostUpdateBegin(
D, INSERT_VALUES, SCATTER_FORWARD);
254 CHKERR VecGhostUpdateEnd(
D, INSERT_VALUES, SCATTER_FORWARD);
268 auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
269 post_proc_fe->generateReferenceElementMesh();
270 post_proc_fe->addFieldValuesPostProc(
"P_REAL");
271 post_proc_fe->addFieldValuesPostProc(
"P_IMAG");
274 CHKERR post_proc_fe->writeFile(
"out_helmholtz.h5m");
284 auto dm = simpleInterface->getDM();
288 CHKERR VecNorm(
D, NORM_2, &nrm2);
290 << std::setprecision(9) <<
"Solution norm " << nrm2;
292 PetscBool test_flg = PETSC_FALSE;
295 constexpr
double regression_test = 97.261672;
296 constexpr
double eps = 1e-6;
297 if (std::abs(nrm2 - regression_test) / regression_test >
eps)
299 "Not converged solution");
305 int main(
int argc,
char *argv[]) {
314 DMType dm_name =
"DMMOFEM";
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
#define CATCH_ERRORS
Catch errors.
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
auto smartCreateDMVector
Get smart vector from DM.
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
boost::ptr_vector< UserDataOperator > & getOpBoundaryRhsPipeline()
Get the Op Boundary Rhs Pipeline object.
#define MOFEM_LOG(channel, severity)
Log.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field.
int main(int argc, char *argv[])
[Check results]
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FTensor::Index< 'k', 3 > k
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)
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
OpCalculateHOJacForFaceImpl< 2 > OpCalculateHOJacForFace
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
DeprecatedCoreInterface Interface
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
const double D
diffusivity
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
Example(MoFEM::Interface &m_field)
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
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.
default operator for EDGE element
friend class UserDataOperator
Set indices on entities on finite element.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MoFEMErrorCode setBoundaryRhsIntegrationRule(RuleHookFun rule)
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator