v0.14.0
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) {}

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 }

◆ 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 }

◆ 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>();
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 
380  CHKERR readMesh();
388 
390 }

◆ 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 }

◆ 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 }

◆ 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 }

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:
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:460
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
WaveEquation::initialCondition
MoFEMErrorCode initialCondition()
Definition: wave_equation.cpp:158
H1
@ H1
continuous field
Definition: definitions.h:85
WaveEquation::setupProblem
MoFEMErrorCode setupProblem()
Definition: wave_equation.cpp:126
OpDomainTimesScalarField
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpDomainTimesScalarField
Definition: wave_equation.cpp:44
OpDomainGradTimesVec
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpGradTimesTensor< 1, 1, SPACE_DIM > OpDomainGradTimesVec
Definition: wave_equation.cpp:46
MoFEM::OpSetHOInvJacToScalarBases< 2 >
Definition: HODataOperators.hpp:78
OpBoundarySource
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< 1, 1 > OpBoundarySource
Definition: wave_equation.cpp:53
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
MoFEM::DMoFEMMeshToLocalVector
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:523
MoFEM::PipelineManager::getBoundaryRhsFE
boost::shared_ptr< FEMethod > & getBoundaryRhsFE()
Definition: PipelineManager.hpp:413
WaveEquation::mField
MoFEM::Interface & mField
Definition: wave_equation.cpp:107
MoFEM::Simple
Simple interface for fast problem set-up.
Definition: Simple.hpp:27
OpBoundaryMass
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpBoundaryMass
Definition: wave_equation.cpp:49
OpDomainMass
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpMass< 1, 1 > OpDomainMass
Definition: wave_equation.cpp:40
order
constexpr int order
Definition: dg_projection.cpp:18
MoFEM::OpCalculateHOJac< 2 >
Definition: HODataOperators.hpp:273
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
WaveEquation::boundaryCondition
MoFEMErrorCode boundaryCondition()
Definition: wave_equation.cpp:163
WaveEquation::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
Definition: wave_equation.cpp:142
OpDomainGradGrad
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< 1, 1, SPACE_DIM > OpDomainGradGrad
Definition: wave_equation.cpp:42
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
simple
void simple(double P1[], double P2[], double P3[], double c[], const int N)
Definition: acoustic.cpp:69
double
wave_speed2
constexpr double wave_speed2
Definition: wave_equation.cpp:55
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
WaveEquation::solveSystem
MoFEMErrorCode solveSystem()
Definition: wave_equation.cpp:272
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
WaveEquation::readMesh
MoFEMErrorCode readMesh()
Definition: wave_equation.cpp:115
WaveEquation::outputResults
MoFEMErrorCode outputResults()
Definition: wave_equation.cpp:369
MoFEM::OpSetHOWeightsOnFace
Modify integration weights on face to take in account higher-order geometry.
Definition: HODataOperators.hpp:122
t
constexpr double t
plate stiffness
Definition: plate.cpp:58
i
FTensor::Index< 'i', SPACE_DIM > i
Definition: hcurl_divergence_operator_2d.cpp:27
MoFEM::OpCalculateScalarFieldValuesDotDot
OpCalculateScalarFieldValuesFromPetscVecImpl< PetscData::CTX_SET_X_TT > OpCalculateScalarFieldValuesDotDot
Definition: UserDataOperators.hpp:275
integration_rule
auto integration_rule
Definition: free_surface.cpp:185
MoFEM::OpUnSetBc
Definition: FormsIntegrators.hpp:49
MoFEM::PipelineManager::getDomainLhsFE
boost::shared_ptr< FEMethod > & getDomainLhsFE()
Definition: PipelineManager.hpp:401
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:308
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition: PetscSmartObj.hpp:221
Monitor
[Push operators to pipeline]
Definition: dynamic_first_order_con_law.cpp:783
approx_order
int approx_order
Definition: test_broken_space.cpp:50
MoFEM::DMMoFEMTSSetMonitor
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:1056
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MoFEM::OpSetBc
Set indices on entities on finite element.
Definition: FormsIntegrators.hpp:38
WaveEquation::boundaryMarker
boost::shared_ptr< std::vector< unsigned char > > boundaryMarker
Definition: wave_equation.cpp:110
MoFEM::OpInvertMatrix
Definition: UserDataOperators.hpp:3249
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:453
OpBoundaryTimeScalarField
FormsIntegrators< BoundaryEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpBaseTimesScalar< 1 > OpBoundaryTimeScalarField
Definition: wave_equation.cpp:51
MoFEM::PetscOptionsGetInt
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
Definition: DeprecatedPetsc.hpp:142
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
WaveEquation::assembleSystem
MoFEMErrorCode assembleSystem()
Definition: wave_equation.cpp:185
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698