v0.13.0
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 /* This file is part of MoFEM.
10  * MoFEM is free software: you can redistribute it and/or modify it under
11  * the terms of the GNU Lesser General Public License as published by the
12  * Free Software Foundation, either version 3 of the License, or (at your
13  * option) any later version.
14  *
15  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
16  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
17  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18  * License for more details.
19  *
20  * You should have received a copy of the GNU Lesser General Public
21  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
22 
23 #include <MoFEM.hpp>
24 
25 using namespace MoFEM;
26 
27 static char help[] = "...\n\n";
28 
29 #include <BasicFiniteElements.hpp>
30 
35 
44 
45 struct Example {
46 
47  Example(MoFEM::Interface &m_field) : mField(m_field) {}
48 
50 
51 private:
52  MoFEM::Interface &mField;
54  boost::shared_ptr<std::vector<unsigned char>> boundaryMarker;
55 
63 
64 };
65 
66 //! [run problem]
69  CHKERR readMesh();
70  CHKERR setupProblem();
71  CHKERR boundaryCondition();
72  CHKERR assembleSystem();
73  CHKERR solveSystem();
74  CHKERR outputResults();
75  CHKERR checkResults();
77 }
78 //! [run problem]
79 
80 //! [Read mesh]
83 
84  CHKERR mField.getInterface(simpleInterface);
85  CHKERR simpleInterface->getOptions();
86  CHKERR simpleInterface->loadFile();
87 
89 }
90 //! [Read mesh]
91 
92 //! [Set up problem]
95  // Add field
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,
104  int order = 6;
105  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
106  CHKERR simpleInterface->setFieldOrder("P_REAL", order);
107  CHKERR simpleInterface->setFieldOrder("P_IMAG", order);
108  CHKERR simpleInterface->setUp();
110 }
111 //! [Set up problem]
112 
113 //! [Applying essential BC]
116 
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);
124  }
125  }
126  // Add vertices to boundary entities
127  Range boundary_vertices;
128  CHKERR mField.get_moab().get_connectivity(boundary_entities,
129  boundary_vertices, true);
130  boundary_entities.merge(boundary_vertices);
131 
132  return boundary_entities;
133  };
134 
135  auto mark_boundary_dofs = [&](Range &&skin_edges) {
136  auto problem_manager = mField.getInterface<ProblemsManager>();
137  auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
138  problem_manager->markDofs(simpleInterface->getProblemName(), ROW,
139  skin_edges, *marker_ptr);
140  return marker_ptr;
141  };
142 
143  auto remove_dofs_from_problem = [&](Range &&skin_edges) {
145  auto problem_manager = mField.getInterface<ProblemsManager>();
146  CHKERR problem_manager->removeDofsOnEntities(
147  simpleInterface->getProblemName(), "P_IMAG", skin_edges, 0, 1);
149  };
150 
151  // Get global local vector of marked DOFs. Is global, since is set for all
152  // DOFs on processor. Is local since only DOFs on processor are in the
153  // vector. To access DOFs use local indices.
154  boundaryMarker = mark_boundary_dofs(get_ents_on_mesh_skin());
155  CHKERR remove_dofs_from_problem(get_ents_on_mesh_skin());
156 
158 }
159 //! [Applying essential BC]
160 
161 //! [Push operators to pipeline]
164  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
165 
166  double k = 90;
167  CHKERR PetscOptionsGetScalar(PETSC_NULL, "", "-k", &k, PETSC_NULL);
168 
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; };
173  auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
174 
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>();
180  pipeline_mng->getOpDomainLhsPipeline().push_back(
181  new OpCalculateHOJacForFace(jac_ptr));
182  pipeline_mng->getOpDomainLhsPipeline().push_back(
183  new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
184  pipeline_mng->getOpDomainLhsPipeline().push_back(
185  new OpSetInvJacH1ForFace(inv_jac_ptr));
186 
187  pipeline_mng->getOpDomainLhsPipeline().push_back(
188  new OpSetBc("P_REAL", true, boundaryMarker));
189 
190  pipeline_mng->getOpDomainLhsPipeline().push_back(
191  new OpDomainGradGrad("P_REAL", "P_REAL", beta));
192  pipeline_mng->getOpDomainLhsPipeline().push_back(
193  new OpDomainGradGrad("P_IMAG", "P_IMAG", beta));
194 
195  pipeline_mng->getOpDomainLhsPipeline().push_back(
196  new OpDomainMass("P_REAL", "P_REAL", k2));
197  pipeline_mng->getOpDomainLhsPipeline().push_back(
198  new OpDomainMass("P_IMAG", "P_IMAG", k2));
199 
200  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
201 
204  };
205 
206  auto set_boundary = [&]() {
208  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
209  new OpSetBc("P_REAL", true, boundaryMarker));
210  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
211  new OpBoundaryMass("P_IMAG", "P_REAL", kp));
212  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
213  new OpBoundaryMass("P_REAL", "P_IMAG", km));
214  pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
215 
216  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
217  new OpSetBc("P_REAL", false, boundaryMarker));
218  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
219  new OpBoundaryMass("P_REAL", "P_REAL", beta));
220  pipeline_mng->getOpBoundaryLhsPipeline().push_back(new OpUnSetBc("P_REAL"));
221 
222  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
223  new OpSetBc("P_REAL", false, boundaryMarker));
224  pipeline_mng->getOpBoundaryRhsPipeline().push_back(
225  new OpBoundarySource("P_REAL", beta));
226  pipeline_mng->getOpBoundaryRhsPipeline().push_back(new OpUnSetBc("P_REAL"));
227 
231  };
232 
233  CHKERR set_domain();
234  CHKERR set_boundary();
235 
237 }
238 //! [Push operators to pipeline]
239 
240 //! [Solve]
243  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
244  auto solver = pipeline_mng->createKSP();
245  CHKERR KSPSetFromOptions(solver);
246  CHKERR KSPSetUp(solver);
247 
248  auto dm = simpleInterface->getDM();
249  auto D = smartCreateDMVector(dm);
250  auto F = smartVectorDuplicate(D);
251 
252  CHKERR KSPSolve(solver, F, D);
253  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
254  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
255  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
257 }
258 //! [Solve]
259 
260 //! [Postprocess results]
263  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
264  pipeline_mng->getDomainLhsFE().reset();
265  pipeline_mng->getDomainRhsFE().reset();
266  pipeline_mng->getBoundaryLhsFE().reset();
267  pipeline_mng->getBoundaryRhsFE().reset();
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");
272  pipeline_mng->getDomainRhsFE() = post_proc_fe;
273  CHKERR pipeline_mng->loopFiniteElements();
274  CHKERR post_proc_fe->writeFile("out_helmholtz.h5m");
276 }
277 //! [Postprocess results]
278 
279 //! [Check results]
282  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
283 
284  auto dm = simpleInterface->getDM();
285  auto D = smartCreateDMVector(dm);
286  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
287  double nrm2;
288  CHKERR VecNorm(D, NORM_2, &nrm2);
289  MOFEM_LOG("WORLD", Sev::inform)
290  << std::setprecision(9) << "Solution norm " << nrm2;
291 
292  PetscBool test_flg = PETSC_FALSE;
293  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-test", &test_flg, PETSC_NULL);
294  if (test_flg) {
295  constexpr double regression_test = 97.261672;
296  constexpr double eps = 1e-6;
297  if (std::abs(nrm2 - regression_test) / regression_test > eps)
298  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
299  "Not converged solution");
300  }
302 }
303 //! [Check results]
304 
305 int main(int argc, char *argv[]) {
306 
307  // Initialisation of MoFEM/PETSc and MOAB data structures
308  const char param_file[] = "param_file.petsc";
309  MoFEM::Core::Initialize(&argc, &argv, param_file, help);
310 
311  try {
312 
313  //! [Register MoFEM discrete manager in PETSc]
314  DMType dm_name = "DMMOFEM";
315  CHKERR DMRegister_MoFEM(dm_name);
316  //! [Register MoFEM discrete manager in PETSc
317 
318  //! [Create MoAB]
319  moab::Core mb_instance; ///< mesh database
320  moab::Interface &moab = mb_instance; ///< mesh database interface
321  //! [Create MoAB]
322 
323  //! [Create MoFEM]
324  MoFEM::Core core(moab); ///< finite element database
325  MoFEM::Interface &m_field = core; ///< finite element database insterface
326  //! [Create MoFEM]
327 
328  //! [Example]
329  Example ex(m_field);
330  CHKERR ex.runProblem();
331  //! [Example]
332  }
333  CATCH_ERRORS;
334 
336 }
std::string param_file
static const double eps
ElementsAndOps< SPACE_DIM >::DomainEle DomainEle
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, FIELD_DIM > OpDomainMass
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_BERNSTEIN_BEZIER_BASE
Definition: definitions.h:77
@ H1
continuous field
Definition: definitions.h:98
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ BLOCKSET
Definition: definitions.h:161
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
auto integration_rule
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:956
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:481
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
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.
Definition: LogManager.hpp:311
#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:305
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
Definition: helmholtz.cpp:39
static char help[]
Definition: helmholtz.cpp:27
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, 2 > OpDomainGradGrad
Definition: helmholtz.cpp:37
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
Definition: helmholtz.cpp:41
FormsIntegrators< EdgeEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: helmholtz.cpp:43
FTensor::Index< 'k', 3 > k
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
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)
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
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]
Definition: plastic.cpp:93
[Example]
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode readMesh()
MoFEMErrorCode checkResults()
MoFEMErrorCode solveSystem()
Example(MoFEM::Interface &m_field)
Definition: helmholtz.cpp:47
Simple * simpleInterface
Definition: helmholtz.cpp:53
MoFEMErrorCode runProblem()
MoFEMErrorCode setupProblem()
MoFEMErrorCode outputResults()
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
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.
Definition: Simple.hpp:33
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator