v0.16.0
Loading...
Searching...
No Matches
Public Member Functions | Private Types | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
Poisson2DHomogeneous Struct Reference
Collaboration diagram for Poisson2DHomogeneous:
[legend]

Public Member Functions

 Poisson2DHomogeneous (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProgram ()
 [Check]
 

Private Types

enum  NORMS { NORM = 0 , LAST_NORM }
 

Private Member Functions

MoFEMErrorCode readMesh ()
 [Read mesh]
 
MoFEMErrorCode setupProblem ()
 [Read mesh]
 
MoFEMErrorCode boundaryCondition ()
 [Setup problem]
 
MoFEMErrorCode assembleSystem ()
 [Boundary condition]
 
MoFEMErrorCode setIntegrationRules ()
 [Assemble system]
 
MoFEMErrorCode solveSystem ()
 [Set integration rules]
 
MoFEMErrorCode outputResults ()
 [Solve system]
 
MoFEMErrorCode checkResults ()
 [Output results]
 

Static Private Member Functions

static double sourceTermFunction (const double x, const double y, const double z)
 

Private Attributes

SimplesimpleInterface
 
MoFEM::InterfacemField
 
int oRder
 
SmartPetscObj< Vec > petscVec
 
int atom_test = 0
 

Detailed Description

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 27 of file poisson_2d_homogeneous.cpp.

Member Enumeration Documentation

◆ NORMS

Constructor & Destructor Documentation

◆ Poisson2DHomogeneous()

Poisson2DHomogeneous::Poisson2DHomogeneous ( MoFEM::Interface m_field)

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode Poisson2DHomogeneous::assembleSystem ( )
private

[Boundary condition]

[Assemble system]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 149 of file poisson_2d_homogeneous.cpp.

149 {
151
152 auto pipeline_mng = mField.getInterface<PipelineManager>();
153
154 { // Push operators to the Pipeline that is responsible for calculating LHS
156 pipeline_mng->getOpDomainLhsPipeline(), {H1});
157 pipeline_mng->getOpDomainLhsPipeline().push_back(
159 }
160
161 { // Push operators to the Pipeline that is responsible for calculating RHS
162
163 auto set_values_to_bc_dofs = [&](auto &fe) {
164 auto get_bc_hook = [&]() {
166 return hook;
167 };
168 fe->preProcessHook = get_bc_hook();
169 };
170
171 // you can skip that if boundary condition is prescribing zero
172 auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
173 using DomainEle =
175 using DomainEleOp = DomainEle::UserDataOperator;
178
179 auto grad_u_vals_ptr = boost::make_shared<MatrixDouble>();
180 pipeline_mng->getOpDomainRhsPipeline().push_back(
182 grad_u_vals_ptr));
183 pipeline_mng->getOpDomainRhsPipeline().push_back(
184 new OpInternal(field_name, grad_u_vals_ptr,
185 [](double, double, double) constexpr { return -1; }));
186 };
187
189 pipeline_mng->getOpDomainRhsPipeline(), {H1});
190 set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
191 calculate_residual_from_set_values_on_bc(
192 pipeline_mng->getOpDomainRhsPipeline());
193 pipeline_mng->getOpDomainRhsPipeline().push_back(
195 }
196
198}
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
#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.
@ PETSC
Standard PETSc assembly.
constexpr auto field_name
constexpr int SPACE_DIM
Add operators pushing bases from local to physical configuration.
Class (Function) to enforce essential constrains.
Definition Essential.hpp:25
Get field gradients at integration pts for scalar field rank 0, i.e. vector field.
Template struct for dimension-specific finite element types.
PipelineManager interface.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
static double sourceTermFunction(const double x, const double y, const double z)

◆ boundaryCondition()

MoFEMErrorCode Poisson2DHomogeneous::boundaryCondition ( )
private

[Setup problem]

[Boundary condition]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 133 of file poisson_2d_homogeneous.cpp.

133 {
135
136 auto bc_mng = mField.getInterface<BcManager>();
137
138 // Remove BCs from blockset name "BOUNDARY_CONDITION" or SETU, note that you
139 // can use regular expression to put list of blocksets;
141 simpleInterface->getProblemName(), "(BOUNDARY_CONDITION|SETU)",
142 std::string(field_name), true);
143
145}
MoFEMErrorCode removeBlockDOFsOnEntities(const std::string problem_name, const std::string block_name, const std::string field_name, int lo, int hi, bool get_low_dim_ents=true, bool is_distributed_mesh=true)
Remove DOFs from problem based on block entities.
Definition BcManager.cpp:72
Boundary condition manager for finite element problem setup.
Template specialization for scalar field boundary conditions.
const std::string getProblemName() const
Get the Problem Name.
Definition Simple.hpp:450

◆ checkResults()

MoFEMErrorCode Poisson2DHomogeneous::checkResults ( )
private

[Output results]

[Check]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 289 of file poisson_2d_homogeneous.cpp.

289 {
291
292 auto check_result_fe_ptr = boost::make_shared<DomainEle>(mField);
293 auto petscVec =
295 (mField.get_comm_rank() == 0) ? LAST_NORM : 0, LAST_NORM);
296
298 check_result_fe_ptr->getOpPtrVector(), {H1})),
299 "Apply transform");
300
301 check_result_fe_ptr->getRuleHook = [](int, int, int p) { return p; };
302 auto analyticalFunction = [&](double x, double y, double z) {
303 return sin(M_PI * x) * sin(M_PI * y);
304 };
305
306 auto u_ptr = boost::make_shared<VectorDouble>();
307
308 check_result_fe_ptr->getOpPtrVector().push_back(
310 auto mValFuncPtr = boost::make_shared<VectorDouble>();
311 check_result_fe_ptr->getOpPtrVector().push_back(
312 new OpGetTensor0fromFunc(mValFuncPtr, analyticalFunction));
313 check_result_fe_ptr->getOpPtrVector().push_back(
314 new OpCalcNormL2Tensor0(u_ptr, petscVec, NORM, mValFuncPtr));
315 CHKERR VecZeroEntries(petscVec);
318 check_result_fe_ptr);
319 CHKERR VecAssemblyBegin(petscVec);
320 CHKERR VecAssemblyEnd(petscVec);
321 MOFEM_LOG_CHANNEL("SELF"); // Clear channel from old tags
322 // print norm in general log
323 if (mField.get_comm_rank() == 0) {
324 const double *norms;
325 CHKERR VecGetArrayRead(petscVec, &norms);
326 MOFEM_TAG_AND_LOG("SELF", Sev::inform, "Errors")
327 << "NORM: " << std::sqrt(norms[NORM]);
328 CHKERR VecRestoreArrayRead(petscVec, &norms);
329 }
330 // compare norm for ctest
331 if (atom_test && !mField.get_comm_rank()) {
332 const double *t_ptr;
333 CHKERR VecGetArrayRead(petscVec, &t_ptr);
334 double ref_norm = 2.2e-04;
335 double cal_norm;
336 switch (atom_test) {
337 case 1: // 2D
338 cal_norm = sqrt(t_ptr[0]);
339 break;
340 default:
341 SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
342 "atom test %d does not exist", atom_test);
343 }
344 if (cal_norm > ref_norm) {
345 SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
346 "atom test %d failed! Calculated Norm %3.16e is greater than "
347 "reference Norm %3.16e",
348 atom_test, cal_norm, ref_norm);
349 }
350 CHKERR VecRestoreArrayRead(petscVec, &t_ptr);
351 }
353}
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
@ MOFEM_ATOM_TEST_INVALID
Definition definitions.h:40
@ MOFEM_INVALID_DATA
Definition definitions.h:36
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Get norm of input VectorDouble for Tensor0.
Specialization for double precision scalar field values calculation.
Get values from scalar function at integration points and save them to VectorDouble for Tensor0.
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition Simple.cpp:799
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition Simple.hpp:429
SmartPetscObj< Vec > petscVec

◆ outputResults()

MoFEMErrorCode Poisson2DHomogeneous::outputResults ( )
private

[Solve system]

[Output results]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 244 of file poisson_2d_homogeneous.cpp.

244 {
246
247 auto pipeline_mng = mField.getInterface<PipelineManager>();
248
249 auto post_proc_fe = boost::make_shared<PostProcFaceEle>(mField);
251 post_proc_fe->getOpPtrVector(), {H1});
252
253 auto u_ptr = boost::make_shared<VectorDouble>();
254 auto grad_u_ptr = boost::make_shared<MatrixDouble>();
255 post_proc_fe->getOpPtrVector().push_back(
257
258 post_proc_fe->getOpPtrVector().push_back(
260
262
263 post_proc_fe->getOpPtrVector().push_back(
264
265 new OpPPMap(post_proc_fe->getPostProcMesh(),
266 post_proc_fe->getMapGaussPts(),
267
268 OpPPMap::DataMapVec{{"U", u_ptr}},
269
270 OpPPMap::DataMapMat{{"GRAD_U", grad_u_ptr}},
271
273
275
276 )
277
278 );
279
280 pipeline_mng->getDomainPostProcFE() = post_proc_fe;
281 CHKERR pipeline_mng->loopFiniteElementsPostProc();
282 CHKERR post_proc_fe->writeFile("out_result.h5m");
283
285}
MoFEMErrorCode loopFiniteElementsPostProc(SmartPetscObj< DM > dm=nullptr)
Iterate postprocessing finite elements.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Post post-proc data at points from hash maps.
std::map< std::string, ScalarDataPtr > DataMapVec
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
boost::shared_ptr< FEMethod > & getDomainPostProcFE()
Get domain postprocessing finite element.

◆ readMesh()

MoFEMErrorCode Poisson2DHomogeneous::readMesh ( )
private

[Read mesh]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 66 of file poisson_2d_homogeneous.cpp.

66 {
68
72 CHKERR mField.getInterface<MeshsetsManager>()->setMeshsetFromFile();
73
75}
Interface for managing meshsets containing materials and boundary conditions.
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition Simple.cpp:191
MoFEMErrorCode getOptions()
get options
Definition Simple.cpp:180

◆ runProgram()

MoFEMErrorCode Poisson2DHomogeneous::runProgram ( )

[Check]

[Run program]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 357 of file poisson_2d_homogeneous.cpp.

357 {
359
368
370}
MoFEMErrorCode setupProblem()
[Read mesh]
MoFEMErrorCode boundaryCondition()
[Setup problem]
MoFEMErrorCode readMesh()
[Read mesh]
MoFEMErrorCode checkResults()
[Output results]
MoFEMErrorCode outputResults()
[Solve system]
MoFEMErrorCode assembleSystem()
[Boundary condition]
MoFEMErrorCode solveSystem()
[Set integration rules]
MoFEMErrorCode setIntegrationRules()
[Assemble system]

◆ setIntegrationRules()

MoFEMErrorCode Poisson2DHomogeneous::setIntegrationRules ( )
private

[Assemble system]

[Set integration rules]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 202 of file poisson_2d_homogeneous.cpp.

202 {
204
205 auto rule_lhs = [](int, int, int p) -> int { return 2 * (p - 1); };
206 auto rule_rhs = [](int, int, int p) -> int { return p; };
207
208 auto pipeline_mng = mField.getInterface<PipelineManager>();
209 CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
210 CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
211
213}
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
Set integration rule for domain left-hand side finite element.

◆ setupProblem()

MoFEMErrorCode Poisson2DHomogeneous::setupProblem ( )
private

[Read mesh]

[Setup problem]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 79 of file poisson_2d_homogeneous.cpp.

79 {
81 Range domain_ents;
82 CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
83 true);
84 auto get_ents_by_dim = [&](const auto dim) {
85 if (dim == SPACE_DIM) {
86 return domain_ents;
87 } else {
88 Range ents;
89 if (dim == 0)
90 CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
91 else
92 CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
93 return ents;
94 }
95 };
96
97 // Select base
98 auto get_base = [&]() {
99 auto domain_ents = get_ents_by_dim(SPACE_DIM);
100 if (domain_ents.empty())
101 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
102 const auto type = type_from_handle(domain_ents[0]);
103 switch (type) {
104 case MBQUAD:
106 case MBHEX:
108 case MBTRI:
110 case MBTET:
112 default:
113 CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type not handled");
114 }
115 return NOBASE;
116 };
117 auto base = get_base();
119
120 int oRder = 3;
121 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-order", &oRder, PETSC_NULLPTR);
122 CHKERR PetscOptionsGetInt(PETSC_NULLPTR, "", "-atom_test", &atom_test,
123 PETSC_NULLPTR);
125
127
129}
std::string type
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition definitions.h:60
@ NOBASE
Definition definitions.h:59
@ DEMKOWICZ_JACOBI_BASE
Definition definitions.h:66
@ H1
continuous field
Definition definitions.h:85
@ MOFEM_NOT_FOUND
Definition definitions.h:33
auto type_from_handle(const EntityHandle h)
get type from entity handle
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
virtual moab::Interface & get_moab()=0
MoFEMErrorCode addDomainField(const std::string name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_ZERO, int verb=-1)
Add field on domain.
Definition Simple.cpp:261
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition Simple.cpp:575
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition Simple.cpp:735

◆ solveSystem()

MoFEMErrorCode Poisson2DHomogeneous::solveSystem ( )
private

[Set integration rules]

[Solve system]

Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 217 of file poisson_2d_homogeneous.cpp.

217 {
219
220 auto pipeline_mng = mField.getInterface<PipelineManager>();
221
222 auto ksp_solver = pipeline_mng->createKSP();
223 CHKERR KSPSetFromOptions(ksp_solver);
224 CHKERR KSPSetUp(ksp_solver);
225
226 // Create RHS and solution vectors
227 auto dm = simpleInterface->getDM();
228 auto F = createDMVector(dm);
229 auto D = vectorDuplicate(F);
230
231 // Solve the system
232 CHKERR KSPSolve(ksp_solver, F, D);
233
234 // Scatter result data on the mesh
235 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
236 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
237 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
238
240}
@ F
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1237
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
double D
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.

◆ sourceTermFunction()

static double Poisson2DHomogeneous::sourceTermFunction ( const double  x,
const double  y,
const double  z 
)
inlinestaticprivate
Examples
mofem/tutorials/scl-1_poisson_2d_homogeneous/poisson_2d_homogeneous.cpp.

Definition at line 52 of file poisson_2d_homogeneous.cpp.

53 {
54 return 2. * M_PI * M_PI * sin(M_PI * x) * sin(M_PI * y);
55 }

Member Data Documentation

◆ atom_test

int Poisson2DHomogeneous::atom_test = 0
private

◆ mField

MoFEM::Interface& Poisson2DHomogeneous::mField
private

◆ oRder

int Poisson2DHomogeneous::oRder
private

◆ petscVec

SmartPetscObj<Vec> Poisson2DHomogeneous::petscVec
private

◆ simpleInterface

Simple* Poisson2DHomogeneous::simpleInterface
private

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