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

Public Member Functions

 WaveEquation (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProgram ()
 

Private Member Functions

MoFEMErrorCode readMesh ()
 
MoFEMErrorCode setupProblem ()
 
MoFEMErrorCode setIntegrationRules ()
 
MoFEMErrorCode initialCondition ()
 
MoFEMErrorCode boundaryCondition ()
 
MoFEMErrorCode assembleSystem ()
 
MoFEMErrorCode solveSystem ()
 
MoFEMErrorCode outputResults ()
 

Private Attributes

MoFEM::InterfacemField
 
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
 

Detailed Description

Examples
wave_equation.cpp.

Definition at line 88 of file wave_equation.cpp.

Constructor & Destructor Documentation

◆ WaveEquation()

WaveEquation::WaveEquation ( MoFEM::Interface m_field)
Examples
wave_equation.cpp.

Definition at line 113 of file wave_equation.cpp.

113: mField(m_field) {}
MoFEM::Interface & mField

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode WaveEquation::assembleSystem ( )
private
Examples
wave_equation.cpp.

Definition at line 185 of file wave_equation.cpp.

185 {
187
188 auto add_domain_base_ops = [&](auto &pipeline) {
189 auto det_ptr = boost::make_shared<VectorDouble>();
190 auto jac_ptr = boost::make_shared<MatrixDouble>();
191 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
192 pipeline.push_back(new OpCalculateHOJac<2>(jac_ptr));
193 pipeline.push_back(new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
194 pipeline.push_back(new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
195 pipeline.push_back(new OpSetHOWeightsOnFace());
196 };
197
198 auto add_domain_lhs_ops = [&](auto &pipeline) {
199 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
200 pipeline.push_back(new OpDomainGradGrad(
201 "U", "U", [](double, double, double) -> double { return 1; }));
202 auto get_c = [this](const double, const double, const double) {
203 auto pipeline_mng = mField.getInterface<PipelineManager>();
204 auto &fe_domain_lhs = pipeline_mng->getDomainLhsFE();
205 return wave_speed2 * fe_domain_lhs->ts_aa;
206 };
207 pipeline.push_back(new OpDomainMass("U", "U", get_c));
208 pipeline.push_back(new OpUnSetBc("U"));
209 };
210
211 auto add_domain_rhs_ops = [&](auto &pipeline) {
212 pipeline.push_back(new OpSetBc("U", true, boundaryMarker));
213 auto grad_u_at_gauss_pts = boost::make_shared<MatrixDouble>();
214 auto dot2_u_at_gauss_pts = boost::make_shared<VectorDouble>();
215 pipeline.push_back(new OpCalculateScalarFieldGradient<SPACE_DIM>(
216 "U", grad_u_at_gauss_pts));
217 pipeline.push_back(
218 new OpCalculateScalarFieldValuesDotDot("U", dot2_u_at_gauss_pts));
219 pipeline.push_back(new OpDomainGradTimesVec(
220 "U", grad_u_at_gauss_pts,
221 [](double, double, double) -> double { return 1; }));
222 pipeline.push_back(new OpDomainTimesScalarField(
223 "U", dot2_u_at_gauss_pts,
224 [](const double, const double, const double) { return wave_speed2; }));
225 pipeline.push_back(new OpUnSetBc("U"));
226 };
227
228 auto add_boundary_base_ops = [&](auto &pipeline) {};
229
230 auto add_lhs_base_ops = [&](auto &pipeline) {
231 pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
232 pipeline.push_back(new OpBoundaryMass(
233 "U", "U", [](const double, const double, const double) { return 1; }));
234 pipeline.push_back(new OpUnSetBc("U"));
235 };
236
237 auto add_rhs_base_ops = [&](auto &pipeline) {
238 pipeline.push_back(new OpSetBc("U", false, boundaryMarker));
239 auto u_at_gauss_pts = boost::make_shared<VectorDouble>();
240 auto boundary_function = [&](const double x, const double y,
241 const double z) {
242 auto pipeline_mng = mField.getInterface<PipelineManager>();
243 auto &fe_rhs = pipeline_mng->getBoundaryRhsFE();
244 const double t = fe_rhs->ts_t;
245 if ((t <= 0.5) && (x < 0.) && (y > -1. / 3) && (y < 1. / 3))
246 return sin(4 * M_PI * t);
247 else
248 return 0.;
249 };
250 pipeline.push_back(new OpCalculateScalarFieldValues("U", u_at_gauss_pts));
251 pipeline.push_back(new OpBoundaryTimeScalarField(
252 "U", u_at_gauss_pts,
253 [](const double, const double, const double) { return 1; }));
254 pipeline.push_back(new OpBoundarySource("U", boundary_function));
255 pipeline.push_back(new OpUnSetBc("U"));
256 };
257
258 auto pipeline_mng = mField.getInterface<PipelineManager>();
259 add_domain_base_ops(pipeline_mng->getOpDomainLhsPipeline());
260 add_domain_base_ops(pipeline_mng->getOpDomainRhsPipeline());
261 add_domain_lhs_ops(pipeline_mng->getOpDomainLhsPipeline());
262 add_domain_rhs_ops(pipeline_mng->getOpDomainRhsPipeline());
263
264 add_boundary_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
265 add_boundary_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
266 add_lhs_base_ops(pipeline_mng->getOpBoundaryLhsPipeline());
267 add_rhs_base_ops(pipeline_mng->getOpBoundaryRhsPipeline());
268
270}
@ 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
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:416
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_TT > OpCalculateScalarFieldValuesDotDot
constexpr double t
plate stiffness
Definition: plate.cpp:59
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Get value at integration points for scalar field.
Set indices on entities on finite element.
Set inverse jacobian to base functions.
Modify integration weights on face to take in account higher-order geometry.
PipelineManager interface.
boost::shared_ptr< FEMethod > & getDomainLhsFE()
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
constexpr double wave_speed2
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass

◆ boundaryCondition()

MoFEMErrorCode WaveEquation::boundaryCondition ( )
private
Examples
wave_equation.cpp.

Definition at line 163 of file wave_equation.cpp.

163 {
165
166 auto bc_mng = mField.getInterface<BcManager>();
167 auto *simple = mField.getInterface<Simple>();
168 CHKERR bc_mng->pushMarkDOFsOnEntities(simple->getProblemName(), "ESSENTIAL",
169 "U", 0, 0);
170
171 auto &bc_map = bc_mng->getBcMapByBlockName();
172 boundaryMarker = boost::make_shared<std::vector<char unsigned>>();
173 for (auto b : bc_map) {
174 if (std::regex_match(b.first, std::regex("(.*)ESSENTIAL(.*)"))) {
175 boundaryMarker->resize(b.second->bcMarkers.size(), 0);
176 for (int i = 0; i != b.second->bcMarkers.size(); ++i) {
177 (*boundaryMarker)[i] |= b.second->bcMarkers[i];
178 }
179 }
180 }
181
183}
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
#define CHKERR
Inline error check.
Definition: definitions.h:535
FTensor::Index< 'i', SPACE_DIM > i
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Simple interface for fast problem set-up.
Definition: Simple.hpp:27

◆ initialCondition()

MoFEMErrorCode WaveEquation::initialCondition ( )
private
Examples
wave_equation.cpp.

Definition at line 158 of file wave_equation.cpp.

158 {
161}

◆ outputResults()

MoFEMErrorCode WaveEquation::outputResults ( )
private
Examples
wave_equation.cpp.

Definition at line 369 of file wave_equation.cpp.

369 {
371
372 // Processes to set output results are integrated in solveSystem()
373
375}

◆ readMesh()

MoFEMErrorCode WaveEquation::readMesh ( )
private
Examples
wave_equation.cpp.

Definition at line 115 of file wave_equation.cpp.

115 {
117
118 auto *simple = mField.getInterface<Simple>();
119 CHKERR mField.getInterface(simple);
120 CHKERR simple->getOptions();
121 CHKERR simple->loadFile();
122
124}

◆ runProgram()

MoFEMErrorCode WaveEquation::runProgram ( )
Examples
wave_equation.cpp.

Definition at line 377 of file wave_equation.cpp.

377 {
379
388
390}
MoFEMErrorCode readMesh()
MoFEMErrorCode setupProblem()
MoFEMErrorCode boundaryCondition()
MoFEMErrorCode initialCondition()
MoFEMErrorCode assembleSystem()
MoFEMErrorCode solveSystem()
MoFEMErrorCode outputResults()
MoFEMErrorCode setIntegrationRules()

◆ setIntegrationRules()

MoFEMErrorCode WaveEquation::setIntegrationRules ( )
private
Examples
wave_equation.cpp.

Definition at line 142 of file wave_equation.cpp.

142 {
144
145 auto integration_rule = [](int o_row, int o_col, int approx_order) {
146 return 2 * approx_order;
147 };
148
149 auto *pipeline_mng = mField.getInterface<PipelineManager>();
150 CHKERR pipeline_mng->setDomainRhsIntegrationRule(integration_rule);
151 CHKERR pipeline_mng->setDomainLhsIntegrationRule(integration_rule);
152 CHKERR pipeline_mng->setBoundaryLhsIntegrationRule(integration_rule);
153 CHKERR pipeline_mng->setBoundaryRhsIntegrationRule(integration_rule);
154
156}
auto integration_rule
static constexpr int approx_order

◆ setupProblem()

MoFEMErrorCode WaveEquation::setupProblem ( )
private
Examples
wave_equation.cpp.

Definition at line 126 of file wave_equation.cpp.

126 {
128
129 auto *simple = mField.getInterface<Simple>();
130 CHKERR simple->addDomainField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
131 CHKERR simple->addBoundaryField("U", H1, AINSWORTH_LEGENDRE_BASE, 1);
132
133 int order = 3;
134 CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &order, PETSC_NULL);
135 CHKERR simple->setFieldOrder("U", order);
136
137 CHKERR simple->setUp();
138
140}
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
constexpr int order
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)

◆ solveSystem()

MoFEMErrorCode WaveEquation::solveSystem ( )
private
Examples
wave_equation.cpp.

Definition at line 272 of file wave_equation.cpp.

272 {
274
275 auto *simple = mField.getInterface<Simple>();
276 auto *pipeline_mng = mField.getInterface<PipelineManager>();
277
278 auto create_post_process_element = [&]() {
279 auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
280
281 auto det_ptr = boost::make_shared<VectorDouble>();
282 auto jac_ptr = boost::make_shared<MatrixDouble>();
283 auto inv_jac_ptr = boost::make_shared<MatrixDouble>();
284 post_proc_fe->getOpPtrVector().push_back(new OpCalculateHOJac<2>(jac_ptr));
285 post_proc_fe->getOpPtrVector().push_back(
286 new OpInvertMatrix<2>(jac_ptr, det_ptr, inv_jac_ptr));
287 post_proc_fe->getOpPtrVector().push_back(
288 new OpSetHOInvJacToScalarBases<2>(H1, inv_jac_ptr));
289
290 auto u_ptr = boost::make_shared<VectorDouble>();
291 post_proc_fe->getOpPtrVector().push_back(
292 new OpCalculateScalarFieldValues("U", u_ptr));
293
295
296 post_proc_fe->getOpPtrVector().push_back(
297
298 new OpPPMap(post_proc_fe->getPostProcMesh(),
299 post_proc_fe->getMapGaussPts(),
300
301 {{"U", u_ptr}},
302
303 {}, {}, {}
304
305 )
306
307 );
308
309 return post_proc_fe;
310 };
311
312 auto set_time_monitor = [&](auto dm, auto solver) {
314 boost::shared_ptr<Monitor> monitor_ptr(
315 new Monitor(dm, create_post_process_element()));
316 boost::shared_ptr<ForcesAndSourcesCore> null;
317 CHKERR DMMoFEMTSSetMonitor(dm, solver, simple->getDomainFEName(),
318 monitor_ptr, null, null);
320 };
321
322 auto set_fieldsplit_preconditioner = [&](auto solver) {
324 SNES snes;
325 CHKERR TSGetSNES(solver, &snes);
326 KSP ksp;
327 CHKERR SNESGetKSP(snes, &ksp);
328 PC pc;
329 CHKERR KSPGetPC(ksp, &pc);
330 PetscBool is_pcfs = PETSC_FALSE;
331 PetscObjectTypeCompare((PetscObject)pc, PCFIELDSPLIT, &is_pcfs);
332
333 if (is_pcfs == PETSC_TRUE) {
334 auto bc_mng = mField.getInterface<BcManager>();
335 auto name_prb = simple->getProblemName();
336 auto is_all_bc = bc_mng->getBlockIS(name_prb, "ESSENTIAL", "U", 0, 0);
337 int is_all_bc_size;
338 CHKERR ISGetSize(is_all_bc, &is_all_bc_size);
339 MOFEM_LOG("EXAMPLE", Sev::inform)
340 << "Field split block size " << is_all_bc_size;
341 CHKERR PCFieldSplitSetIS(pc, PETSC_NULL,
342 is_all_bc); // boundary block
343 }
345 };
346
347 auto dm = simple->getDM();
348 auto D = createDMVector(dm);
349 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_FORWARD);
350
351 auto ts = pipeline_mng->createTSIM2();
352
353 CHKERR TSSetSolution(ts, D);
354 auto DD = vectorDuplicate(D);
355 CHKERR TS2SetSolution(ts, D, DD);
356 CHKERR set_time_monitor(dm, ts);
357 CHKERR TSSetFromOptions(ts);
358 CHKERR set_fieldsplit_preconditioner(ts);
359 CHKERR TSSetUp(ts);
360 CHKERR TSSolve(ts, NULL);
361
362 CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
363 CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
364 CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
365
367}
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:447
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:440
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode)
set local (or ghosted) vector values on mesh for partition only
Definition: DMMoFEM.cpp:509
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1003
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
double D
PetscErrorCode DMMoFEMTSSetMonitor(DM dm, TS ts, const std::string fe_name, boost::shared_ptr< MoFEM::FEMethod > method, boost::shared_ptr< MoFEM::BasicMethod > pre_only, boost::shared_ptr< MoFEM::BasicMethod > post_only)
Set Monitor To TS solver.
Definition: DMMoFEM.cpp:1042
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Post post-proc data at points from hash maps.
[Push operators to pipeline]

Member Data Documentation

◆ boundaryMarker

boost::shared_ptr<std::vector<unsigned char> > WaveEquation::boundaryMarker
private
Examples
wave_equation.cpp.

Definition at line 110 of file wave_equation.cpp.

◆ mField

MoFEM::Interface& WaveEquation::mField
private
Examples
wave_equation.cpp.

Definition at line 107 of file wave_equation.cpp.


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