v0.9.2
lesson5_helmholtz.cpp
Go to the documentation of this file.
1 /**
2  * \file lesson5_helmholtz.cpp
3  * \example lesson5_helmholtz.cpp
4  *
5  * Using PipelineManager interface calculate the divergence of base functions, and
6  * integral of flux on the boundary. Since the h-div space is used, volume
7  * integral and boundary integral should give the same result.
8  */
9 
10 /* This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
23 
24 #include <MoFEM.hpp>
25 
26 using namespace MoFEM;
27 
28 static char help[] = "...\n\n";
29 
30 #include <BasicFiniteElements.hpp>
31 
33 
38 
39 #include <BaseOps.hpp>
40 
45 
46 struct Example {
47 
48  Example(MoFEM::Interface &m_field) : mField(m_field) {}
49 
50  MoFEMErrorCode runProblem();
51 
52 private:
53  MoFEM::Interface &mField;
54 
55  MoFEMErrorCode setUP();
56  MoFEMErrorCode bC();
57  MoFEMErrorCode OPs();
58  MoFEMErrorCode kspSolve();
59  MoFEMErrorCode postProcess();
60 
62 };
63 
66  CHKERR setUP();
67  CHKERR bC();
68  CHKERR OPs();
69  CHKERR kspSolve();
70  CHKERR postProcess();
72 }
73 
74 //! [Set up problem]
77  Simple *simple = mField.getInterface<Simple>();
78  // Add field
79  CHKERR simple->addDomainField("U_REAL", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
80  1);
81  CHKERR simple->addDomainField("U_IMAG", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
82  1);
83  CHKERR simple->addBoundaryField("U_REAL", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
84  1);
85  CHKERR simple->addBoundaryField("U_IMAG", H1, AINSWORTH_BERNSTEIN_BEZIER_BASE,
86  1);
87  constexpr int order = 10;
88  CHKERR simple->setFieldOrder("U_REAL", order);
89  CHKERR simple->setFieldOrder("U_IMAG", order);
90  CHKERR simple->setUp();
92 }
93 //! [Set up problem]
94 
95 //! [Applying essential BC]
99 }
100 //! [Applying essential BC]
101 
102 //! [Push operators to pipeline]
105  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
106 
107  auto vol_source_function = [](const double x, const double y,
108  const double z) {
109  const double xc = 0;
110  const double yc = -1.25;
111  const double xs = x - xc;
112  const double ys = y - yc;
113  constexpr double eps = 60;
114  return exp(-pow(eps * sqrt(xs * xs + ys * ys), 2));
115  };
116 
117  constexpr double k = 90;
118 
119  auto beta = [](const double, const double, const double) { return -1; };
120  auto k2 = [k](const double, const double, const double) { return pow(k, 2); };
121  auto kp = [k](const double, const double, const double) { return k; };
122  auto km = [k](const double, const double, const double) { return -k; };
123  auto integration_rule = [](int, int, int p_data) { return 2 * p_data; };
124 
125  auto set_domain = [&]() {
127  pipeline_mng->getOpDomainLhsPipeline().push_back(
129  pipeline_mng->getOpDomainLhsPipeline().push_back(new OpSetInvJacH1ForFace(invJac));
130  pipeline_mng->getOpDomainLhsPipeline().push_back(
131  new OpDomainGradGrad("U_REAL", "U_REAL", beta));
132  pipeline_mng->getOpDomainLhsPipeline().push_back(
133  new OpDomainGradGrad("U_IMAG", "U_IMAG", beta));
134 
135  pipeline_mng->getOpDomainLhsPipeline().push_back(
136  new OpDomainMass("U_REAL", "U_REAL", k2));
137  pipeline_mng->getOpDomainLhsPipeline().push_back(
138  new OpDomainMass("U_IMAG", "U_IMAG", k2));
139 
140  pipeline_mng->getOpDomainRhsPipeline().push_back(
141  new OpDomainSource("U_REAL", vol_source_function));
142 
143  CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
144  CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
146  };
147 
148  auto set_boundary = [&]() {
150  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
151  new OpBoundaryMass("U_IMAG", "U_REAL", kp));
152  pipeline_mng->getOpBoundaryLhsPipeline().push_back(
153  new OpBoundaryMass("U_REAL", "U_IMAG", km));
154  CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
156  };
157 
158  CHKERR set_domain();
159  CHKERR set_boundary();
160 
162 }
163 //! [Push operators to pipeline]
164 
165 //! [Solve]
168  Simple *simple = mField.getInterface<Simple>();
169  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
170  auto solver = pipeline_mng->createKSP();
171  CHKERR KSPSetFromOptions(solver);
172  CHKERR KSPSetUp(solver);
173 
174  auto dm = simple->getDM();
175  auto D = smartCreateDMVector(dm);
176  auto F = smartVectorDuplicate(D);
177 
178  CHKERR KSPSolve(solver, F, D);
179  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
180  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
181  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
183 }
184 //! [Solve]
185 
186 //! [Postprocess results]
189  PipelineManager *pipeline_mng = mField.getInterface<PipelineManager>();
190  pipeline_mng->getDomainLhsFE().reset();
191  pipeline_mng->getDomainRhsFE().reset();
192  pipeline_mng->getBoundaryLhsFE().reset();
193  pipeline_mng->getBoundaryRhsFE().reset();
194  auto post_proc_fe = boost::make_shared<PostProcFaceOnRefinedMesh>(mField);
195  post_proc_fe->generateReferenceElementMesh();
196  post_proc_fe->addFieldValuesPostProc("U_REAL");
197  post_proc_fe->addFieldValuesPostProc("U_IMAG");
198  pipeline_mng->getDomainRhsFE() = post_proc_fe;
199  CHKERR pipeline_mng->loopFiniteElements();
200  CHKERR post_proc_fe->writeFile("out_helmholtz.h5m");
202 }
203 //! [Postprocess results]
204 
205 int main(int argc, char *argv[]) {
206 
207  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
208 
209  try {
210 
211  //! [Register MoFEM discrete manager in PETSc]
212  DMType dm_name = "DMMOFEM";
213  CHKERR DMRegister_MoFEM(dm_name);
214  //! [Register MoFEM discrete manager in PETSc
215 
216  //! [Create MoAB]
217  moab::Core mb_instance; ///< mesh database
218  moab::Interface &moab = mb_instance; ///< mesh database interface
219  //! [Create MoAB]
220 
221  //! [Create MoFEM]
222  MoFEM::Core core(moab); ///< finite element database
223  MoFEM::Interface &m_field = core; ///< finite element database insterface
224  //! [Create MoFEM]
225 
226  //! [Load mesh]
227  Simple *simple = m_field.getInterface<Simple>();
228  CHKERR simple->getOptions();
229  CHKERR simple->loadFile();
230  //! [Load mesh]
231 
232  //! [Example]
233  Example ex(m_field);
234  CHKERR ex.runProblem();
235  //! [Example]
236  }
237  CATCH_ERRORS;
238 
240 }
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::ptr_vector< UserDataOperator > & getOpDomainRhsPipeline()
Get the Op Domain Rhs Pipeline object.
Deprecated interface functions.
MoFEMErrorCode kspSolve()
[Push operators to pipeline]
OpTools< DomainEleOp >::OpSource< 2 > OpDomainSource
[Grad residual operator]
Definition: BaseOps.hpp:239
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:445
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:75
MoFEMErrorCode setBoundaryLhsIntegrationRule(RuleHookFun rule)
MatrixDouble invJac
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:483
Example(MoFEM::Interface &m_field)
OpTools< DomainEleOp >::OpMass OpDomainMass
Core (interface) class.
Definition: Core.hpp:50
DataForcesAndSourcesCore::EntData EntData
SmartPetscObj< Vec > smartVectorDuplicate(SmartPetscObj< Vec > &vec)
Create duplicate vector of smart vector.
Definition: AuxPETSc.hpp:238
static char help[]
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
boost::ptr_vector< UserDataOperator > & getOpDomainLhsPipeline()
Get the Op Domain Lhs Pipeline object.
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:51
Edge finite elementUser is implementing own operator at Gauss points level, by own object derived fro...
auto smartCreateDMVector
Get smart vector from DM.
Definition: DMMoFEM.hpp:939
OpTools< EdgeEleOp >::OpMass OpBoundaryMass
OpTools< DomainEleOp >::OpGradGrad< 2 > OpDomainGradGrad
MoFEMErrorCode runProblem()
[Operator]
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
MoFEMErrorCode OPs()
[Boundary condition]
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Calculate inverse of jacobian for face element.
MoFEMErrorCode bC()
[Create common data]
[Source operator]
Definition: BaseOps.hpp:97
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
boost::shared_ptr< FEMethod > & getDomainRhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
PipelineManager interface.
#define CHKERR
Inline error check.
Definition: definitions.h:602
boost::shared_ptr< FEMethod > & getBoundaryLhsFE()
MoFEMErrorCode setUP()
[Set up problem]
MoFEMErrorCode loopFiniteElements(SmartPetscObj< DM > dm=nullptr)
Iterate finite elements.
Transform local reference derivatives of shape functions to global derivatives.
constexpr int order
FTensor::Index< 'k', 2 > k
Definition: ContactOps.hpp:28
MoFEMErrorCode postProcess()
[Integrate]
MoFEMErrorCode setDomainRhsIntegrationRule(RuleHookFun rule)
int main(int argc, char *argv[])
[Postprocess results]
boost::ptr_vector< UserDataOperator > & getOpBoundaryLhsPipeline()
Get the Op Boundary Lhs Pipeline object.
[Source operator]
Definition: BaseOps.hpp:137
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1879
Data on single entity (This is passed as argument to DataOperator::doWork)
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:413
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:439
continuous field
Definition: definitions.h:177
static const double eps
MoFEMErrorCode setDomainLhsIntegrationRule(RuleHookFun rule)
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:69
Face finite elementUser is implementing own operator at Gauss point level, by own object derived from...