v0.13.1
helmholtz.cpp
Go to the documentation of this file.
1/**
2 * \file helmholtz.cpp
3 * \example helmholtz.cpp
4 *
5 * Using PipelineManager interface calculate Helmholtz problem. Example show how
6 * to manage complex variable fields.
7 */
8
9
10
11#include <MoFEM.hpp>
12
13using namespace MoFEM;
14
15static char help[] = "...\n\n";
16
17#include <BasicFiniteElements.hpp>
18
23
24constexpr int SPACE_DIM = 2;
25
34
35
36
37struct Example {
38
39 Example(MoFEM::Interface &m_field) : mField(m_field) {}
40
42
43private:
46 boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
47
55
56};
57
58//! [run problem]
69}
70//! [run problem]
71
72//! [Read mesh]
75
79
81}
82//! [Read mesh]
83
84//! [Set up problem]
87 // Add field
96 int order = 6;
97 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
102}
103//! [Set up problem]
104
105//! [Applying essential BC]
108
109 auto get_ents_on_mesh_skin = [&]() {
110 Range boundary_entities;
112 std::string entity_name = it->getName();
113 if (entity_name.compare(0, 2, "BC") == 0) {
114 CHKERR it->getMeshsetIdEntitiesByDimension(mField.get_moab(), 1,
115 boundary_entities, true);
116 }
117 }
118 // Add vertices to boundary entities
119 Range boundary_vertices;
120 CHKERR mField.get_moab().get_connectivity(boundary_entities,
121 boundary_vertices, true);
122 boundary_entities.merge(boundary_vertices);
123
124 return boundary_entities;
125 };
126
127 auto mark_boundary_dofs = [&](Range &&skin_edges) {
128 auto problem_manager = mField.getInterface<ProblemsManager>();
129 auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
130 problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
131 skin_edges, *marker_ptr);
132 return marker_ptr;
133 };
134
135 auto remove_dofs_from_problem = [&](Range &&skin_edges) {
137 auto problem_manager = mField.getInterface<ProblemsManager>();
138 CHKERR problem_manager->removeDofsOnEntities(
139 simpleInterface->getProblemName(), "P_IMAG", skin_edges, 0, 1);
141 };
142
143 // Get global local vector of marked DOFs. Is global, since is set for all
144 // DOFs on processor. Is local since only DOFs on processor are in the
145 // vector. To access DOFs use local indices.
146 boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
147 CHKERR remove_dofs_from_problem(get_ents_on_mesh_skin());
148
150}
151//! [Applying essential BC]
152
153//! [Push operators to pipeline]
157
158 double k = 90;
159 CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-k", &k, PETSC_NULL);
160
161 auto beta = [](const double, const double, const double) { return -1; };
162 auto k2 = [k](const double, const double, const double) { return pow(k, 2); };
163 auto kp = [k](const double, const double, const double) { return k; };
164 auto km = [k](const double, const double, const double) { return -k; };
165 auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
166
167 auto set_domain = [&]() {
169 auto det_ptr = boost::make_shared<VectorDouble>();
170 auto jac_ptr = boost::make_shared<MatrixDouble>();
171 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
172 pipeline_mng->getOpDomainLhsPipeline().push_back(
173 new OpCalculateHOJac<2>(jac_ptr));
174 pipeline_mng->getOpDomainLhsPipeline().push_back(
175 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
176 pipeline_mng->getOpDomainLhsPipeline().push_back(
177 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
178 pipeline_mng->getOpDomainLhsPipeline().push_back(
180
181 pipeline_mng->getOpDomainLhsPipeline().push_back(
182 new OpSetBc("P_REAL", true, boundaryMarker));
183
184 pipeline_mng->getOpDomainLhsPipeline().push_back(
185 new OpDomainGradGrad("P_REAL", "P_REAL", beta));
186 pipeline_mng->getOpDomainLhsPipeline().push_back(
187 new OpDomainGradGrad("P_IMAG", "P_IMAG", beta));
188
189 pipeline_mng->getOpDomainLhsPipeline().push_back(
190 new OpDomainMass("P_REAL", "P_REAL", k2));
191 pipeline_mng->getOpDomainLhsPipeline().push_back(
192 new OpDomainMass("P_IMAG", "P_IMAG", k2));
193
194 pipeline_mng->getOpDomainLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
195
198 };
199
200 auto set_boundary = [&]() {
202 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
203 new OpSetBc("P_REAL", true, boundaryMarker));
204 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
205 new OpBoundaryMass("P_IMAG", "P_REAL", kp));
206 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
207 new OpBoundaryMass("P_REAL", "P_IMAG", km));
208 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
209
210 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
211 new OpSetBc("P_REAL", false, boundaryMarker));
212 pipeline_mng->getOpBoundaryLhsPipeline().push_back(
213 new OpBoundaryMass("P_REAL", "P_REAL", beta));
214 pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
215
216 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
217 new OpSetBc("P_REAL", false, boundaryMarker));
218 pipeline_mng->getOpBoundaryRhsPipeline().push_back(
219 new OpBoundarySource("P_REAL", beta));
220 pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("P_REAL"));
221
225 };
226
227 CHKERR set_domain();
228 CHKERR set_boundary();
229
231}
232//! [Push operators to pipeline]
233
234//! [Solve]
238 auto solver = pipeline_mng->createKSP();
239 CHKERR KSPSetFromOptions(solver);
240 CHKERR KSPSetUp(solver);
241
242 auto dm = simpleInterface->getDM();
243 auto D = smartCreateDMVector(dm);
244 auto F = smartVectorDuplicate(D);
245
246 CHKERR KSPSolve(solver, F, D);
247 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
248 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
249 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
251}
252//! [Solve]
253
254//! [Postprocess results]
258 pipeline_mng->getDomainLhsFE().reset();
259 pipeline_mng->getDomainRhsFE().reset();
260 pipeline_mng->getBoundaryLhsFE().reset();
261 pipeline_mng->getBoundaryRhsFE().reset();
262
264
265 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
266
267 auto p_real_ptr = boost::make_shared<VectorDouble>();
268 auto p_imag_ptr = boost::make_shared<VectorDouble>();
269
270 post_proc_fe->getOpPtrVector().push_back(
271 new OpCalculateScalarFieldValues("P_REAL", p_real_ptr));
272 post_proc_fe->getOpPtrVector().push_back(
273 new OpCalculateScalarFieldValues("P_IMAG", p_imag_ptr));
274
276
277 post_proc_fe->getOpPtrVector().push_back(
278
279 new OpPPMap(
280
281 post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
282
283 {{"P_REAL", p_real_ptr}, {"P_IMAG", p_imag_ptr}},
284
285 {}, {}, {}
286
287 )
288
289 );
290
291 pipeline_mng->getDomainRhsFE() = post_proc_fe;
292 CHKERR pipeline_mng->loopFiniteElements();
293 CHKERR post_proc_fe->writeFile("out_helmholtz.h5m");
295}
296//! [Postprocess results]
297
298//! [Check results]
302
303 auto dm = simpleInterface->getDM();
304 auto D = smartCreateDMVector(dm);
305 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
306 double nrm2;
307 CHKERR VecNorm(D, NORM_2, &nrm2);
308 MOFEM_LOG("WORLD", Sev::inform)
309 << std::setprecision(9) << "Solution norm " << nrm2;
310
311 PetscBool test_flg = PETSC_FALSE;
312 CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
313 if (test_flg) {
314 constexpr double regression_test = 97.261672;
315 constexpr double eps = 1e-6;
316 if (std::abs(nrm2 - regression_test) / regression_test > eps)
317 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
318 "Not converged solution");
319 }
321}
322//! [Check results]
323
324int main(int argc, char *argv[]) {
325
326 // Initialisation of MoFEM/PETSc and MOAB data structures
327 const char param_file[] = "param_file.petsc";
329
330 try {
331
332 //! [Register MoFEM discrete manager in PETSc]
333 DMType dm_name = "DMMOFEM";
334 CHKERR DMRegister_MoFEM(dm_name);
335 //! [Register MoFEM discrete manager in PETSc
336
337 //! [Create MoAB]
338 moab::Core mb_instance; ///< mesh database
339 moab::Interface &moab = mb_instance; ///< mesh database interface
340 //! [Create MoAB]
341
342 //! [Create MoFEM]
343 MoFEM::Core core(moab); ///< finite element database
344 MoFEM::Interface &m_field = core; ///< finite element database insterface
345 //! [Create MoFEM]
346
347 //! [Example]
348 Example ex(m_field);
349 CHKERR ex.runProblem();
350 //! [Example]
351 }
353
355}
std::string param_file
ForcesAndSourcesCore::UserDataOperator UserDataOperator
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
@ ROW
Definition: definitions.h:123
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:372
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:64
@ H1
continuous field
Definition: definitions.h:85
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:346
@ BLOCKSET
Definition: definitions.h:148
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
#define CHKERR
Inline error check.
Definition: definitions.h:535
auto integration_rule
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMMoFEM.cpp:470
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:47
auto smartCreateDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:965
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.
Definition: LogManager.hpp:301
#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]
Definition: helmholtz.cpp:324
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
Definition: helmholtz.cpp:29
static char help[]
Definition: helmholtz.cpp:15
constexpr int SPACE_DIM
Definition: helmholtz.cpp:24
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: helmholtz.cpp:27
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
Definition: helmholtz.cpp:31
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:33
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:56
implementation of Data Operators for Forces and Sources
Definition: MoFEM.hpp:24
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)
PetscErrorCode PetscOptionsGetScalar(PetscOptions *, const char pre[], const char name[], PetscScalar *dval, PetscBool *set)
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
CoreTmp< 0 > Core
Definition: Core.hpp:1086
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1955
const double D
diffusivity
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, SPACE_DIM > OpBoundaryMass
[Only used with Hencky/nonlinear material]
Definition: plastic.cpp:68
[Example]
Definition: plastic.cpp:120
MoFEMErrorCode boundaryCondition()
[Set up problem]
Definition: helmholtz.cpp:106
MoFEMErrorCode assembleSystem()
[Applying essential BC]
Definition: helmholtz.cpp:154
MoFEMErrorCode readMesh()
[run problem]
Definition: helmholtz.cpp:73
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
[Push operators to pipeline]
Definition: helmholtz.cpp:235
Example(MoFEM::Interface &m_field)
Definition: helmholtz.cpp:39
Simple * simpleInterface
Definition: helmholtz.cpp:45
MoFEMErrorCode runProblem()
MoFEM::Interface & mField
Definition: plastic.cpp:127
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: plastic.cpp:146
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
[Solve]
Definition: helmholtz.cpp:255
virtual moab::Interface & get_moab()=0
Core (interface) class.
Definition: Core.hpp:82
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:72
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:112
Deprecated interface functions.
Get value at integration points for scalar field.
Post post-proc data at points from hash maps.
Set indices on entities on finite element.
Modify integration weights on face to take in account higher-order geometry.
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.
Definition: Simple.hpp:26
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:374
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:290
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:806
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:304
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:589
MoFEMErrorCode addBoundaryField(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 boundary.
Definition: Simple.cpp:392
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:749
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:334
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.