v0.14.0
Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
Electrostatics Struct Reference
Collaboration diagram for Electrostatics:
[legend]

Public Member Functions

 Electrostatics (MoFEM::Interface &m_field)
 
MoFEMErrorCode runProgram ()
 [Get Charges] More...
 

Private Types

enum  VecElements { ZERO = 0, ONE = 1, LAST_ELEMENT }
 

Private Member Functions

MoFEMErrorCode readMesh ()
 [Read mesh] More...
 
MoFEMErrorCode setupProblem ()
 [Read mesh] More...
 
MoFEMErrorCode boundaryCondition ()
 [Setup problem] More...
 
MoFEMErrorCode setIntegrationRules ()
 [Boundary condition] More...
 
MoFEMErrorCode assembleSystem ()
 [Set integration rules] More...
 
MoFEMErrorCode solveSystem ()
 [Assemble system] More...
 
MoFEMErrorCode outputResults ()
 [Solve system] More...
 
MoFEMErrorCode getTotalEnergy ()
 [Output results] More...
 
MoFEMErrorCode getElectrodeCharge ()
 [Get Total Energy] More...
 

Private Attributes

int oRder
 
MoFEM::InterfacemField
 
SimplesimpleInterface
 
std::string domainField
 
boost::shared_ptr< std::map< int, BlockData > > permBlockSetsPtr
 
boost::shared_ptr< std::map< int, BlockData > > intBlockSetsPtr
 
boost::shared_ptr< std::map< int, BlockData > > electrodeBlockSetsPtr
 
boost::shared_ptr< DataAtIntegrationPtscommonDataPtr
 
boost::shared_ptr< ForcesAndSourcesCoreinterFaceRhsFe
 
boost::shared_ptr< ForcesAndSourcesCoreelectrodeRhsFe
 
double aLpha = 0.0
 
double bEta = 0.0
 
SmartPetscObj< Vec > petscVec
 
SmartPetscObj< Vec > petscVecEnergy
 
PetscBool out_skin = PETSC_FALSE
 
PetscBool is_partitioned = PETSC_FALSE
 
int atomTest = 0
 

Detailed Description

Definition at line 12 of file electrostatics.cpp.

Member Enumeration Documentation

◆ VecElements

Enumerator
ZERO 
ONE 
LAST_ELEMENT 

Definition at line 50 of file electrostatics.cpp.

50 { ZERO = 0, ONE = 1, LAST_ELEMENT };

Constructor & Destructor Documentation

◆ Electrostatics()

Electrostatics::Electrostatics ( MoFEM::Interface m_field)

Definition at line 54 of file electrostatics.cpp.

55  : domainField("POTENTIAL"), mField(m_field) {}

Member Function Documentation

◆ assembleSystem()

MoFEMErrorCode Electrostatics::assembleSystem ( )
private

[Set integration rules]

[Assemble system]

Definition at line 324 of file electrostatics.cpp.

324  {
326 
327  auto pipeline_mng = mField.getInterface<PipelineManager>();
328 
329  commonDataPtr = boost::make_shared<DataAtIntegrationPts>(mField);
330  auto add_base_ops = [&](auto &pipeline) {
332 
333  pipeline.push_back(
335  };
336 
337  add_base_ops(pipeline_mng->getOpDomainLhsPipeline());
338  auto epsilon = [&](const double, const double, const double) {
339  return commonDataPtr->blockPermittivity;
340  };
341 
342  { // Push operators to the Pipeline that is responsible for calculating LHS
343  pipeline_mng->getOpDomainLhsPipeline().push_back(
345  }
346 
347  { // Push operators to the Pipeline that is responsible for calculating RHS
348  auto set_values_to_bc_dofs = [&](auto &fe) {
349  auto get_bc_hook = [&]() {
351  return hook;
352  };
353  fe->preProcessHook = get_bc_hook();
354  };
355  // Set essential BC
356  auto calculate_residual_from_set_values_on_bc = [&](auto &pipeline) {
357  using OpInternal =
360 
361  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
362  pipeline_mng->getOpDomainRhsPipeline().push_back(
364  grad_u_ptr));
365 
366  add_base_ops(pipeline_mng->getOpDomainRhsPipeline());
367 
368  auto minus_epsilon = [&](double, double, double) constexpr {
369  return -commonDataPtr->blockPermittivity;
370  };
371  pipeline_mng->getOpDomainRhsPipeline().push_back(
372  new OpInternal(domainField, grad_u_ptr, minus_epsilon));
373  };
374 
375  set_values_to_bc_dofs(pipeline_mng->getDomainRhsFE());
376  calculate_residual_from_set_values_on_bc(
377  pipeline_mng->getOpDomainRhsPipeline());
378  auto bodySourceTerm = [&](const double, const double, const double) {
379  return bodySource;
380  };
381  pipeline_mng->getOpDomainRhsPipeline().push_back(
382  new OpBodySourceVectorb(domainField, bodySourceTerm));
383  }
384 
385  interFaceRhsFe = boost::shared_ptr<ForcesAndSourcesCore>(
387  interFaceRhsFe->getRuleHook = [](int, int, int p) { return p; };
388 
389  {
390  interFaceRhsFe->getOpPtrVector().push_back(
392 
393  auto sIgma = [&](const double, const double, const double) {
394  return commonDataPtr->blockChrgDens;
395  };
396 
397  interFaceRhsFe->getOpPtrVector().push_back(
398  new OpInterfaceRhsVectorF(domainField, sIgma));
399  }
400 
402 }

◆ boundaryCondition()

MoFEMErrorCode Electrostatics::boundaryCondition ( )
private

[Setup problem]

[Boundary condition]

Definition at line 294 of file electrostatics.cpp.

294  {
296 
297  auto bc_mng = mField.getInterface<BcManager>();
298 
299  // Remove_BCs_from_blockset name "BOUNDARY_CONDITION";
300  CHKERR bc_mng->removeBlockDOFsOnEntities<BcScalarMeshsetType<BLOCKSET>>(
301  simpleInterface->getProblemName(), "BOUNDARY_CONDITION",
302  std::string(domainField), true);
303 
305 }

◆ getElectrodeCharge()

MoFEMErrorCode Electrostatics::getElectrodeCharge ( )
private

[Get Total Energy]

[Get Charges]

Definition at line 606 of file electrostatics.cpp.

606  {
608  auto op_loop_side = new OpLoopSide<SideEle>(
610 
612  op_loop_side->getOpPtrVector(), {H1});
613 
614  auto grad_u_ptr_charge = boost::make_shared<MatrixDouble>();
615  auto e_ptr_charge = boost::make_shared<MatrixDouble>();
616 
617  op_loop_side->getOpPtrVector().push_back(
619  grad_u_ptr_charge));
620 
621  op_loop_side->getOpPtrVector().push_back(
622  new OpElectricField(e_ptr_charge, grad_u_ptr_charge));
623  auto d_jump = boost::make_shared<MatrixDouble>();
624  op_loop_side->getOpPtrVector().push_back(new OpElectricDispJump<SPACE_DIM>(
625  domainField, e_ptr_charge, d_jump, commonDataPtr, permBlockSetsPtr));
626 
627  electrodeRhsFe = boost::shared_ptr<ForcesAndSourcesCore>(
629  electrodeRhsFe->getRuleHook = [](int, int, int p) { return p; };
630 
631  // push all the operators in on the side to the electrodeRhsFe
632  electrodeRhsFe->getOpPtrVector().push_back(op_loop_side);
633 
634  electrodeRhsFe->getOpPtrVector().push_back(new OpElectrodeCharge<SPACE_DIM>(
636  CHKERR VecZeroEntries(petscVec);
638  simpleInterface->getDM(), "ELECTRODE", electrodeRhsFe,
640  CHKERR VecAssemblyBegin(petscVec);
641  CHKERR VecAssemblyEnd(petscVec);
642 
643  if (!mField.get_comm_rank()) {
644  const double *array;
645 
646  CHKERR(VecGetArrayRead(petscVec, &array));
647  double aLpha = array[0]; // Use explicit index instead of ZERO
648  double bEta = array[1]; // Use explicit index instead of ONE
649  MOFEM_LOG_CHANNEL("SELF");
650  MOFEM_LOG_C("SELF", Sev::inform,
651  "CHARGE_ELEC_1: %6.15f , CHARGE_ELEC_2: %6.15f", aLpha, bEta);
652 
653  CHKERR(VecRestoreArrayRead(petscVec, &array));
654  }
655 
656  if (atomTest && !mField.get_comm_rank()) {
657  double cal_charge_elec1;
658  double cal_charge_elec2;
659  double cal_total_energy;
660  const double *c_ptr, *te_ptr;
661 
662  // Get a pointer to the PETSc vector data
663  CHKERR(VecGetArrayRead(petscVec, &c_ptr));
664  CHKERR(VecGetArrayRead(petscVecEnergy, &te_ptr));
665 
666  // Expected charges at the electrodes
667  double ref_charge_elec1 = 50.0;
668  double ref_charge_elec2 = -50.0;
669  // Expected total energy of the system
670  double ref_tot_energy = 500.0;
671 
672  switch (atomTest) {
673  case 1: // 2D & 3D test
674  cal_charge_elec1 = c_ptr[0]; // Read the charge at the first electrode
675  cal_charge_elec2 = c_ptr[1]; // Read the charge at the second electrode
676  cal_total_energy = te_ptr[0]; // Read the total energy of the system
677  break;
678  default:
679  SETERRQ1(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
680  "atom test %d does not exist", atomTest);
681  }
682 
683  // Validate the results
684  if (std::abs(ref_charge_elec1 - cal_charge_elec1) > 1e-10 ||
685  std::abs(ref_charge_elec2 - cal_charge_elec2) > 1e-10 ||
686  std::abs(ref_tot_energy - cal_total_energy) > 1e-10) {
687  SETERRQ1(
688  PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
689  "atom test %d failed! Calculated values do not match expected values",
690  atomTest);
691  }
692 
693  CHKERR(VecRestoreArrayRead(petscVec,
694  &c_ptr)); // Restore the PETSc vector array
695  CHKERR(VecRestoreArrayRead(petscVecEnergy, &te_ptr));
696  }
697 
699 }

◆ getTotalEnergy()

MoFEMErrorCode Electrostatics::getTotalEnergy ( )
private

[Output results]

[Get Total Energy]

Definition at line 539 of file electrostatics.cpp.

539  {
541  auto pip_energy = mField.getInterface<PipelineManager>();
542  pip_energy->getDomainRhsFE().reset();
543  pip_energy->getDomainLhsFE().reset();
544  pip_energy->getBoundaryLhsFE().reset();
545  pip_energy->getBoundaryRhsFE().reset();
546  pip_energy->getOpDomainRhsPipeline().clear();
547  pip_energy->getOpDomainLhsPipeline().clear();
548 
549  // gets the map of the internal domain entity range to get the total energy
550 
551  boost::shared_ptr<std::map<int, BlockData>> intrnlDomnBlckSetPtr =
552  boost::make_shared<std::map<int, BlockData>>();
553  Range internal_domain; // range of entities marked the internal domain
555  if (bit->getName().compare(0, 10, "DOMAIN_INT") == 0) {
556  const int id = bit->getMeshsetId();
557  auto &block_data = (*intrnlDomnBlckSetPtr)[id];
558 
559  CHKERR mField.get_moab().get_entities_by_dimension(
560  bit->getMeshset(), SPACE_DIM, block_data.internalDomainEnts, true);
561  internal_domain.merge(block_data.internalDomainEnts);
562  block_data.iD = id;
563  }
564  }
565  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
566  internal_domain);
567 
569  pip_energy->getOpDomainRhsPipeline(), {H1});
570 
571 
572  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
573  auto e_field_ptr = boost::make_shared<MatrixDouble>();
574 
575  pip_energy->getOpDomainRhsPipeline().push_back(
577  pip_energy->getOpDomainRhsPipeline().push_back(
578  new OpElectricField(e_field_ptr, grad_u_ptr));
579 
580  commonDataPtr = boost::make_shared<DataAtIntegrationPts>(mField);
581 
582  pip_energy->getOpDomainRhsPipeline().push_back(
583  new OpTotalEnergy(domainField, grad_u_ptr, permBlockSetsPtr,
584  intrnlDomnBlckSetPtr, commonDataPtr, petscVecEnergy));
585 
586  pip_energy->loopFiniteElements();
587  CHKERR VecAssemblyBegin(petscVecEnergy);
588  CHKERR VecAssemblyEnd(petscVecEnergy);
589 
590  double total_energy = 0.0; // declaration for total energy
591  if (!mField.get_comm_rank()) {
592  const double *array;
593 
594  CHKERR VecGetArrayRead(petscVecEnergy, &array);
595  total_energy = array[ZERO];
596  MOFEM_LOG_CHANNEL("SELF");
597  MOFEM_LOG_C("SELF", Sev::inform, "Total Energy: %6.15f", total_energy);
598  CHKERR VecRestoreArrayRead(petscVecEnergy, &array);
599  }
600 
602 }

◆ outputResults()

MoFEMErrorCode Electrostatics::outputResults ( )
private

[Solve system]

[Output results]

Definition at line 449 of file electrostatics.cpp.

449  {
451  auto pipeline_mng = mField.getInterface<PipelineManager>();
452  pipeline_mng->getDomainRhsFE().reset();
453  pipeline_mng->getBoundaryLhsFE().reset();
454  pipeline_mng->getBoundaryRhsFE().reset(); //
455  pipeline_mng->getDomainLhsFE().reset();
456  auto post_proc_fe = boost::make_shared<PostProcEle>(mField);
457  auto u_ptr = boost::make_shared<VectorDouble>();
458  auto grad_u_ptr = boost::make_shared<MatrixDouble>();
459  auto e_field_ptr = boost::make_shared<MatrixDouble>();
460  // lamda function to calculate electric field
461  auto add_efield_ops = [&](auto &pipeline) {
462  // add higher order operator
464  // calculate field values
465  pipeline.push_back(new OpCalculateScalarFieldValues(domainField, u_ptr));
466  // calculate gradient
467  pipeline.push_back(
469  // calculate electric field
470  pipeline.push_back(new OpElectricField(e_field_ptr, grad_u_ptr));
471  };
472 
473  add_efield_ops(post_proc_fe->getOpPtrVector());
474 
475  auto e_field_times_perm_ptr = boost::make_shared<MatrixDouble>();
476  auto energy_density_ptr = boost::make_shared<VectorDouble>();
477 
478  post_proc_fe->getOpPtrVector().push_back(
479  new OpGradTimesPerm(domainField, e_field_ptr, e_field_times_perm_ptr,
481  post_proc_fe->getOpPtrVector().push_back(
482  new OpEnergyDensity(domainField, e_field_ptr, energy_density_ptr,
484 
486  post_proc_fe->getOpPtrVector().push_back(new OpPPMap(
487  post_proc_fe->getPostProcMesh(), post_proc_fe->getMapGaussPts(),
488 
489  OpPPMap::DataMapVec{{"POTENTIAL", u_ptr},
490  {"ENERGY_DENSITY", energy_density_ptr}},
492  {"ELECTRIC_FIELD", e_field_ptr},
493  {"ELECTRIC_DISPLACEMENT", e_field_times_perm_ptr},
494  },
496 
498 
499  );
500 
501  pipeline_mng->getDomainRhsFE() = post_proc_fe;
502  CHKERR pipeline_mng->loopFiniteElements();
503  CHKERR post_proc_fe->writeFile("out.h5m");
504 
505  if (out_skin && SPACE_DIM == 3) {
506  auto post_proc_skin = boost::make_shared<PostProcFaceEle>(mField);
507  auto op_loop_skin = new OpLoopSide<SideEle>(
509 
510  add_efield_ops(op_loop_skin->getOpPtrVector());
511 
512  op_loop_skin->getOpPtrVector().push_back(
513  new OpGradTimesPerm(domainField, e_field_ptr, e_field_times_perm_ptr,
515  op_loop_skin->getOpPtrVector().push_back(
516  new OpEnergyDensity(domainField, e_field_ptr, energy_density_ptr,
518 
519  // push op to boundary element
520  post_proc_skin->getOpPtrVector().push_back(op_loop_skin);
521 
522  post_proc_skin->getOpPtrVector().push_back(new OpPPMap(
523  post_proc_skin->getPostProcMesh(), post_proc_skin->getMapGaussPts(),
524  OpPPMap::DataMapVec{{"POTENTIAL", u_ptr},
525  {"ENERGY_DENSITY", energy_density_ptr}},
526  OpPPMap::DataMapMat{{"ELECTRIC_FIELD", e_field_ptr},
527  {"ELECTRIC_DISPLACEMENT", e_field_times_perm_ptr}},
529 
531  post_proc_skin);
532  CHKERR post_proc_skin->writeFile("out_skin.h5m");
533  }
535 }

◆ readMesh()

MoFEMErrorCode Electrostatics::readMesh ( )
private

[Read mesh]

Definition at line 58 of file electrostatics.cpp.

58  {
62 
64  true; // create lower dimensional element to ensure sharing of partitioed
65  // entities at the boundaries.
68 }

◆ runProgram()

MoFEMErrorCode Electrostatics::runProgram ( )

[Get Charges]

[Run program]

Definition at line 703 of file electrostatics.cpp.

◆ setIntegrationRules()

MoFEMErrorCode Electrostatics::setIntegrationRules ( )
private

[Boundary condition]

[Set integration rules]

Definition at line 309 of file electrostatics.cpp.

309  {
311 
312  auto rule_lhs = [](int, int, int p) -> int { return 2 * p; };
313  auto rule_rhs = [](int, int, int p) -> int { return p; };
314 
315  auto pipeline_mng = mField.getInterface<PipelineManager>();
316  CHKERR pipeline_mng->setDomainLhsIntegrationRule(rule_lhs);
317  CHKERR pipeline_mng->setDomainRhsIntegrationRule(rule_rhs);
318 
320 }

◆ setupProblem()

MoFEMErrorCode Electrostatics::setupProblem ( )
private

[Read mesh]

[Setup problem]

Definition at line 72 of file electrostatics.cpp.

72  {
74 
75  Range domain_ents;
76  CHKERR mField.get_moab().get_entities_by_dimension(0, SPACE_DIM, domain_ents,
77  true);
78  auto get_ents_by_dim = [&](const auto dim) {
79  if (dim == SPACE_DIM) {
80  return domain_ents;
81  } else {
82  Range ents;
83  if (dim == 0)
84  CHKERR mField.get_moab().get_connectivity(domain_ents, ents, true);
85  else
86  CHKERR mField.get_moab().get_entities_by_dimension(0, dim, ents, true);
87  return ents;
88  }
89  };
90 
91  // Select base for the field based on the element type
92  auto get_base = [&]() {
93  auto domain_ents = get_ents_by_dim(SPACE_DIM);
94  if (domain_ents.empty())
95  CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Empty mesh");
96  const auto type = type_from_handle(domain_ents[0]);
97  switch (type) {
98  case MBQUAD:
99  return DEMKOWICZ_JACOBI_BASE;
100  case MBHEX:
101  return DEMKOWICZ_JACOBI_BASE;
102  case MBTRI:
104  case MBTET:
106  default:
107  CHK_THROW_MESSAGE(MOFEM_NOT_FOUND, "Element type is not handled");
108  }
109  return NOBASE;
110  };
111 
112  auto base = get_base();
115 
116  int oRder = 2; // default order
117 
118  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-order", &oRder, PETSC_NULL);
120  CHKERR PetscOptionsGetInt(PETSC_NULL, "", "-atomTest", &atomTest, PETSC_NULL);
121  commonDataPtr = boost::make_shared<DataAtIntegrationPts>(mField);
122 
123  // gets the map of the permittivity attributes and the block sets
124  permBlockSetsPtr = boost::make_shared<std::map<int, BlockData>>();
125  Range mat_electr_ents; // range of entities with the permittivity
127  if (bit->getName().compare(0, 12, "MAT_ELECTRIC") == 0) {
128  const int id = bit->getMeshsetId();
129  auto &block_data = (*permBlockSetsPtr)[id];
130 
131  CHKERR mField.get_moab().get_entities_by_dimension(
132  bit->getMeshset(), SPACE_DIM, block_data.domainEnts, true);
133  mat_electr_ents.merge(block_data.domainEnts);
134 
135  std::vector<double> attributes;
136  bit->getAttributes(attributes);
137  if (attributes.size() < 1) {
138  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
139  " At least one permittivity attributes should be given but "
140  "found %d",
141  attributes.size());
142  }
143  block_data.iD = id; // id of the block
144  block_data.epsPermit = attributes[0]; // permittivity value of the block
145  }
146  }
147 
148  // gets the map of the charge attributes and the block sets
149  intBlockSetsPtr = boost::make_shared<std::map<int, BlockData>>();
150  Range int_electr_ents; // range of entities with the charge
152  if (bit->getName().compare(0, 12, "INT_ELECTRIC") == 0) {
153  const int id = bit->getMeshsetId();
154  auto &block_data = (*intBlockSetsPtr)[id];
155 
156  CHKERR mField.get_moab().get_entities_by_dimension(
157  bit->getMeshset(), SPACE_DIM - 1, block_data.interfaceEnts, true);
158  int_electr_ents.merge(block_data.interfaceEnts);
159 
160  std::vector<double> attributes;
161  bit->getAttributes(attributes);
162  if (attributes.size() < 1) {
163  SETERRQ1(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
164  "At least one charge attributes should be given but found %d",
165  attributes.size());
166  }
167 
168  block_data.iD = id; // id of the block
169  block_data.chargeDensity = attributes[0]; // charge value of the block
170  }
171  }
172  // gets the map of the electrode entity range in the block sets
173  electrodeBlockSetsPtr = boost::make_shared<std::map<int, BlockData>>();
174  Range electrode_ents; // range of entities with the electrode
175  int electrodeCount = 0;
177  if (bit->getName().compare(0, 9, "ELECTRODE") == 0) {
178  const int id = bit->getMeshsetId();
179  auto &block_data = (*electrodeBlockSetsPtr)[id];
180  ++electrodeCount;
181 
182  CHKERR mField.get_moab().get_entities_by_dimension(
183  bit->getMeshset(), SPACE_DIM - 1, block_data.electrodeEnts, true);
184  electrode_ents.merge(block_data.electrodeEnts);
185  block_data.iD = id;
186  auto print_range_on_procs = [&](const std::string &name, int meshsetId,
187  const Range &range) {
188  MOFEM_LOG("SYNC", Sev::inform)
189  << name << " in meshID: " << id << " with range "
190  << block_data.electrodeEnts << " on proc ["
191  << mField.get_comm_rank() << "] \n";
193  };
194  print_range_on_procs(bit->getName(), id, electrode_ents);
195  if (electrodeCount > 2) {
196  SETERRQ(PETSC_COMM_SELF, MOFEM_INVALID_DATA,
197  "Three or more electrode blocksets found");
198  ;
199  }
200  }
201  }
202 
203  // sync entities
204  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
205  mat_electr_ents);
206  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
207  int_electr_ents);
208  CHKERR mField.getInterface<CommInterface>()->synchroniseEntities(
209  electrode_ents);
210 
211  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-out_skin", &out_skin,
212  PETSC_NULL);
213  CHKERR PetscOptionsGetBool(PETSC_NULL, "", "-is_partitioned", &is_partitioned,
214  PETSC_NULL);
215  // get the skin entities
216  Skinner skinner(&mField.get_moab());
217  Range skin_tris;
218  CHKERR skinner.find_skin(0, mat_electr_ents, false, skin_tris);
219  Range proc_skin;
220  ParallelComm *pcomm =
221  ParallelComm::get_pcomm(&mField.get_moab(), MYPCOMM_INDEX);
222  if (is_partitioned) {
223  CHKERR pcomm->filter_pstatus(skin_tris,
224  PSTATUS_SHARED | PSTATUS_MULTISHARED,
225  PSTATUS_NOT, -1, &proc_skin);
226  } else {
227  proc_skin = skin_tris;
228  }
229  // add the skin entities to the field
232  "SKIN");
236  // add the interface entities to the field
237  CHKERR mField.add_finite_element("INTERFACE");
242  SPACE_DIM - 1, "INTERFACE");
243  // add the electrode entities to the field
244  CHKERR mField.add_finite_element("ELECTRODE");
249  "ELECTRODE");
250 
251  // sync field entities
252  mField.getInterface<CommInterface>()->synchroniseFieldEntities(domainField);
254  CHKERR simpleInterface->defineProblem(PETSC_TRUE);
259  CHKERR mField.build_finite_elements("INTERFACE");
260  CHKERR mField.build_finite_elements("ELECTRODE");
261 
265 
266  DMType dm_name = "DMMOFEM";
267  CHKERR DMRegister_MoFEM(dm_name);
268 
270  dm = createDM(mField.get_comm(), dm_name);
271 
272  // create dm instance
273  CHKERR DMSetType(dm, dm_name);
274 
276 
277  // initialise petsc vector for required processor
278  int local_size;
279  if (mField.get_comm_rank() == 0) // get_comm_rank() gets processor number
280 
281  local_size = LAST_ELEMENT; // last element gives size of vector
282 
283  else
284  // other processors (e.g. 1, 2, 3, etc.)
285  local_size = 0; // local size of vector is zero on other processors
286 
290 }

◆ solveSystem()

MoFEMErrorCode Electrostatics::solveSystem ( )
private

[Assemble system]

[Solve system]

< Null element does

Definition at line 406 of file electrostatics.cpp.

406  {
408 
409  auto pipeline_mng = mField.getInterface<PipelineManager>();
410 
411  auto ksp_solver = pipeline_mng->createKSP();
412 
413  boost::shared_ptr<ForcesAndSourcesCore> null; ///< Null element does
415  interFaceRhsFe, null, null);
416 
417  CHKERR KSPSetFromOptions(ksp_solver);
418  CHKERR KSPSetUp(ksp_solver);
419 
420  // Create RHS and solution vectors
421  auto dm = simpleInterface->getDM();
422  auto F = createDMVector(dm);
423  auto D = vectorDuplicate(F);
424  // Solve the system
425  CHKERR KSPSolve(ksp_solver, F, D);
426 
427  CHKERR VecGhostUpdateBegin(F, INSERT_VALUES, SCATTER_FORWARD);
428  CHKERR VecGhostUpdateEnd(F, INSERT_VALUES, SCATTER_FORWARD);
429 
430  double fnorm;
431  CHKERR VecNorm(F, NORM_2, &fnorm);
432  CHKERR PetscPrintf(PETSC_COMM_WORLD, "F norm = %9.8e\n", fnorm);
433 
434  // Scatter result data on the mesh
435  CHKERR VecGhostUpdateBegin(D, INSERT_VALUES, SCATTER_FORWARD);
436  CHKERR VecGhostUpdateEnd(D, INSERT_VALUES, SCATTER_FORWARD);
437 
438  double dnorm;
439  CHKERR VecNorm(D, NORM_2, &dnorm);
440  CHKERR PetscPrintf(PETSC_COMM_WORLD, "D norm = %9.8e\n", dnorm);
441 
442  CHKERR DMoFEMMeshToLocalVector(dm, D, INSERT_VALUES, SCATTER_REVERSE);
443 
445 }

Member Data Documentation

◆ aLpha

double Electrostatics::aLpha = 0.0
private

Definition at line 44 of file electrostatics.cpp.

◆ atomTest

int Electrostatics::atomTest = 0
private

Definition at line 51 of file electrostatics.cpp.

◆ bEta

double Electrostatics::bEta = 0.0
private

Definition at line 45 of file electrostatics.cpp.

◆ commonDataPtr

boost::shared_ptr<DataAtIntegrationPts> Electrostatics::commonDataPtr
private

Definition at line 38 of file electrostatics.cpp.

◆ domainField

std::string Electrostatics::domainField
private

Definition at line 34 of file electrostatics.cpp.

◆ electrodeBlockSetsPtr

boost::shared_ptr<std::map<int, BlockData> > Electrostatics::electrodeBlockSetsPtr
private

Definition at line 37 of file electrostatics.cpp.

◆ electrodeRhsFe

boost::shared_ptr<ForcesAndSourcesCore> Electrostatics::electrodeRhsFe
private

Definition at line 41 of file electrostatics.cpp.

◆ intBlockSetsPtr

boost::shared_ptr<std::map<int, BlockData> > Electrostatics::intBlockSetsPtr
private

Definition at line 36 of file electrostatics.cpp.

◆ interFaceRhsFe

boost::shared_ptr<ForcesAndSourcesCore> Electrostatics::interFaceRhsFe
private

Definition at line 40 of file electrostatics.cpp.

◆ is_partitioned

PetscBool Electrostatics::is_partitioned = PETSC_FALSE
private

Definition at line 49 of file electrostatics.cpp.

◆ mField

MoFEM::Interface& Electrostatics::mField
private

Definition at line 32 of file electrostatics.cpp.

◆ oRder

int Electrostatics::oRder
private

Definition at line 31 of file electrostatics.cpp.

◆ out_skin

PetscBool Electrostatics::out_skin = PETSC_FALSE
private

Definition at line 48 of file electrostatics.cpp.

◆ permBlockSetsPtr

boost::shared_ptr<std::map<int, BlockData> > Electrostatics::permBlockSetsPtr
private

Definition at line 35 of file electrostatics.cpp.

◆ petscVec

SmartPetscObj<Vec> Electrostatics::petscVec
private

Definition at line 46 of file electrostatics.cpp.

◆ petscVecEnergy

SmartPetscObj<Vec> Electrostatics::petscVecEnergy
private

Definition at line 47 of file electrostatics.cpp.

◆ simpleInterface

Simple* Electrostatics::simpleInterface
private

Definition at line 33 of file electrostatics.cpp.


The documentation for this struct was generated from the following file:
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition: UnknownInterface.hpp:93
Electrostatics::readMesh
MoFEMErrorCode readMesh()
[Read mesh]
Definition: electrostatics.cpp:58
MoFEM::PipelineManager::getDomainRhsFE
boost::shared_ptr< FEMethod > & getDomainRhsFE()
Definition: PipelineManager.hpp:405
MYPCOMM_INDEX
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
H1
@ H1
continuous field
Definition: definitions.h:85
MoFEM::DMMoFEMKSPSetComputeRHS
PetscErrorCode DMMoFEMKSPSetComputeRHS(DM dm, const char fe_name[], MoFEM::FEMethod *method, MoFEM::BasicMethod *pre_only, MoFEM::BasicMethod *post_only)
set KSP right hand side evaluation function
Definition: DMMoFEM.cpp:637
MOFEM_LOG_CHANNEL
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:284
Electrostatics::atomTest
int atomTest
Definition: electrostatics.cpp:51
MoFEM::OpPostProcMapInMoab::DataMapVec
std::map< std::string, boost::shared_ptr< VectorDouble > > DataMapVec
Definition: PostProcBrokenMeshInMoabBase.hpp:700
OpBodySourceVectorb
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpBodySourceVectorb
Definition: electrostatics.hpp:45
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::modify_finite_element_add_field_row
virtual MoFEMErrorCode modify_finite_element_add_field_row(const std::string &fe_name, const std::string name_row)=0
set field row which finite element use
NOBASE
@ NOBASE
Definition: definitions.h:59
CHK_THROW_MESSAGE
#define CHK_THROW_MESSAGE(err, msg)
Check and throw MoFEM exception.
Definition: definitions.h:609
IntElementForcesAndSourcesCore
intPostProc< SPACE_DIM >::intEle IntElementForcesAndSourcesCore
Definition: electrostatics.hpp:35
Electrostatics::out_skin
PetscBool out_skin
Definition: electrostatics.cpp:48
MoFEM::Simple::buildProblem
MoFEMErrorCode buildProblem()
Build problem.
Definition: Simple.cpp:669
OpInterfaceRhsVectorF
FormsIntegrators< IntEleOp >::Assembly< PETSC >::LinearForm< GAUSS >::OpSource< BASE_DIM, FIELD_DIM > OpInterfaceRhsVectorF
Definition: electrostatics.hpp:43
OpElectrodeCharge
Definition: electrostatics.hpp:357
MoFEM::CoreInterface::get_comm_rank
virtual int get_comm_rank() const =0
MoFEM::Simple::loadFile
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name, LoadFileFunc loadFunc=defaultLoadFileFunc)
Load mesh file.
Definition: Simple.cpp:194
MoFEM::PETSC
@ PETSC
Definition: FormsIntegrators.hpp:105
MoFEM::PipelineManager
PipelineManager interface.
Definition: PipelineManager.hpp:24
Electrostatics::commonDataPtr
boost::shared_ptr< DataAtIntegrationPts > commonDataPtr
Definition: electrostatics.cpp:38
Electrostatics::permBlockSetsPtr
boost::shared_ptr< std::map< int, BlockData > > permBlockSetsPtr
Definition: electrostatics.cpp:35
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
Electrostatics::outputResults
MoFEMErrorCode outputResults()
[Solve system]
Definition: electrostatics.cpp:449
Electrostatics::bEta
double bEta
Definition: electrostatics.cpp:45
Electrostatics::petscVecEnergy
SmartPetscObj< Vec > petscVecEnergy
Definition: electrostatics.cpp:47
MoFEM::DMMoFEMAddElement
PetscErrorCode DMMoFEMAddElement(DM dm, std::string fe_name)
add element to dm
Definition: DMMoFEM.cpp:497
Electrostatics::getTotalEnergy
MoFEMErrorCode getTotalEnergy()
[Output results]
Definition: electrostatics.cpp:539
MoFEM::OpCalculateScalarFieldGradient
Get field gradients at integration pts for scalar filed rank 0, i.e. vector field.
Definition: UserDataOperators.hpp:1293
MoFEM::Simple::getOptions
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:180
OpElectricDispJump
Definition: electrostatics.hpp:300
MoFEM::Simple::getAddSkeletonFE
bool & getAddSkeletonFE()
Get the addSkeletonFE.
Definition: Simple.hpp:452
MoFEM::Simple::getDM
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:747
CHKERR
#define CHKERR
Inline error check.
Definition: definitions.h:548
Electrostatics::mField
MoFEM::Interface & mField
Definition: electrostatics.cpp:32
MoFEM::createDMVector
auto createDMVector(DM dm)
Get smart vector from DM.
Definition: DMMoFEM.hpp:1099
MoFEM::CoreInterface::add_finite_element
virtual MoFEMErrorCode add_finite_element(const std::string &fe_name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
add finite element
bodySource
const double bodySource
Definition: electrostatics.hpp:11
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
OpTotalEnergy
Definition: electrostatics.hpp:135
Electrostatics::domainField
std::string domainField
Definition: electrostatics.cpp:34
MoFEM::DMoFEMLoopFiniteElementsUpAndLowRank
PetscErrorCode DMoFEMLoopFiniteElementsUpAndLowRank(DM dm, const char fe_name[], MoFEM::FEMethod *method, int low_rank, int up_rank, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:567
MoFEM::CoreInterface::modify_finite_element_add_field_col
virtual MoFEMErrorCode modify_finite_element_add_field_col(const std::string &fe_name, const std::string name_row)=0
set field col which finite element use
Electrostatics::boundaryCondition
MoFEMErrorCode boundaryCondition()
[Setup problem]
Definition: electrostatics.cpp:294
SPACE_DIM
constexpr int SPACE_DIM
Definition: child_and_parent.cpp:16
Electrostatics::electrodeRhsFe
boost::shared_ptr< ForcesAndSourcesCore > electrodeRhsFe
Definition: electrostatics.cpp:41
MoFEM::createDM
auto createDM(MPI_Comm comm, const std::string dm_type_name)
Creates smart DM object.
Definition: PetscSmartObj.hpp:141
OpBlockChargeDensity
Definition: electrostatics.hpp:74
MoFEM::BcManager
Simple interface for fast problem set-up.
Definition: BcManager.hpp:25
Electrostatics::oRder
int oRder
Definition: electrostatics.cpp:31
MoFEM::CoreInterface::build_finite_elements
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elements.
MOFEM_LOG_C
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:311
double
bit
auto bit
set bit
Definition: hanging_node_approx.cpp:75
convert.type
type
Definition: convert.py:64
MoFEM::FormsIntegrators::Assembly
Assembly methods.
Definition: FormsIntegrators.hpp:317
OpPPMap
OpPostProcMapInMoab< SPACE_DIM, SPACE_DIM > OpPPMap
Definition: photon_diffusion.cpp:29
MoFEM::BcScalarMeshsetType
Definition: BcManager.hpp:15
OpDomainLhsMatrixK
FormsIntegrators< DomainEleOp >::Assembly< PETSC >::BiLinearForm< GAUSS >::OpGradGrad< BASE_DIM, FIELD_DIM, SPACE_DIM > OpDomainLhsMatrixK
Definition: electrostatics.hpp:40
MOFEM_LOG_SYNCHRONISE
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:345
MoFEM::OpCalculateScalarFieldValues
Get value at integration points for scalar field.
Definition: UserDataOperators.hpp:82
MoFEM::Simple::addDomainField
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:264
MoFEM::PipelineManager::createKSP
SmartPetscObj< KSP > createKSP(SmartPetscObj< DM > dm=nullptr)
Create KSP (linear) solver.
Definition: PipelineManager.cpp:49
OpGradTimesPerm
Definition: electrostatics.hpp:247
MoFEM::DMRegister_MoFEM
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMoFEM.cpp:43
OpElectricField
Definition: electrostatics.hpp:406
MoFEM::GAUSS
@ GAUSS
Definition: FormsIntegrators.hpp:136
MoFEM::Simple::getDomainFEName
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:368
MoFEM::type_from_handle
auto type_from_handle(const EntityHandle h)
get type from entity handle
Definition: Templates.hpp:1898
Electrostatics::intBlockSetsPtr
boost::shared_ptr< std::map< int, BlockData > > intBlockSetsPtr
Definition: electrostatics.cpp:36
MoFEM::CoreInterface::add_ents_to_finite_element_by_dim
virtual MoFEMErrorCode add_ents_to_finite_element_by_dim(const EntityHandle entities, const int dim, const std::string &name, const bool recursive=true)=0
add entities to finite element
OpEnergyDensity
Definition: electrostatics.hpp:195
MoFEM::CoreInterface::modify_finite_element_add_field_data
virtual MoFEMErrorCode modify_finite_element_add_field_data(const std::string &fe_name, const std::string name_field)=0
set finite element field data
MoFEM::OpPostProcMapInMoab::DataMapMat
std::map< std::string, boost::shared_ptr< MatrixDouble > > DataMapMat
Definition: PostProcBrokenMeshInMoabBase.hpp:701
MOFEM_NOT_FOUND
@ MOFEM_NOT_FOUND
Definition: definitions.h:33
OpGradTimesTensor
FormsIntegrators< DomainEleOp >::Assembly< A >::LinearForm< I >::OpGradTimesTensor< 1, FIELD_DIM, SPACE_DIM > OpGradTimesTensor
Definition: operators_tests.cpp:48
Electrostatics::is_partitioned
PetscBool is_partitioned
Definition: electrostatics.cpp:49
MoFEM::Simple::setFieldOrder
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:545
MoFEM::Simple::buildFiniteElements
MoFEMErrorCode buildFiniteElements()
Build finite elements.
Definition: Simple.cpp:643
MoFEM::AddHOOps
Add operators pushing bases from local to physical configuration.
Definition: HODataOperators.hpp:413
MoFEM::EssentialPreProc< TemperatureCubitBcData >
Specialization for TemperatureCubitBcData.
Definition: EssentialTemperatureCubitBcData.hpp:91
Range
MF_ZERO
@ MF_ZERO
Definition: definitions.h:111
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
DEMKOWICZ_JACOBI_BASE
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:66
MoFEM::CommInterface
Managing BitRefLevels.
Definition: CommInterface.hpp:21
_IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_
#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.
Definition: MeshsetsManager.hpp:71
Electrostatics::setIntegrationRules
MoFEMErrorCode setIntegrationRules()
[Boundary condition]
Definition: electrostatics.cpp:309
Electrostatics::simpleInterface
Simple * simpleInterface
Definition: electrostatics.cpp:33
Electrostatics::solveSystem
MoFEMErrorCode solveSystem()
[Assemble system]
Definition: electrostatics.cpp:406
BLOCKSET
@ BLOCKSET
Definition: definitions.h:161
Electrostatics::setupProblem
MoFEMErrorCode setupProblem()
[Read mesh]
Definition: electrostatics.cpp:72
OpBlockPermittivity
Definition: electrostatics.hpp:103
Electrostatics::petscVec
SmartPetscObj< Vec > petscVec
Definition: electrostatics.cpp:46
AINSWORTH_LEGENDRE_BASE
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:60
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:31
MoFEM::Simple::defineFiniteElements
MoFEMErrorCode defineFiniteElements()
Define finite elements.
Definition: Simple.cpp:452
Electrostatics::ONE
@ ONE
Definition: electrostatics.cpp:50
MoFEM::createVectorMPI
auto createVectorMPI(MPI_Comm comm, PetscInt n, PetscInt N)
Create MPI Vector.
Definition: PetscSmartObj.hpp:202
MoFEM::Simple::buildFields
MoFEMErrorCode buildFields()
Build fields.
Definition: Simple.cpp:554
Electrostatics::interFaceRhsFe
boost::shared_ptr< ForcesAndSourcesCore > interFaceRhsFe
Definition: electrostatics.cpp:40
Electrostatics::getElectrodeCharge
MoFEMErrorCode getElectrodeCharge()
[Get Total Energy]
Definition: electrostatics.cpp:606
MoFEM::Simple::addBoundaryField
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:354
ReactionDiffusionEquation::D
const double D
diffusivity
Definition: reaction_diffusion.cpp:20
Electrostatics::assembleSystem
MoFEMErrorCode assembleSystem()
[Set integration rules]
Definition: electrostatics.cpp:324
MOFEM_ATOM_TEST_INVALID
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:40
MoFEM::Simple::defineProblem
MoFEMErrorCode defineProblem(const PetscBool is_partitioned=PETSC_TRUE)
define problem
Definition: Simple.cpp:524
MoFEM::SmartPetscObj< DM >
Electrostatics::electrodeBlockSetsPtr
boost::shared_ptr< std::map< int, BlockData > > electrodeBlockSetsPtr
Definition: electrostatics.cpp:37
Electrostatics::ZERO
@ ZERO
Definition: electrostatics.cpp:50
MoFEM::Simple::getProblemName
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:389
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMoFEM.cpp:586
convert.int
int
Definition: convert.py:64
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
Electrostatics::LAST_ELEMENT
@ LAST_ELEMENT
Definition: electrostatics.cpp:50
Electrostatics::aLpha
double aLpha
Definition: electrostatics.cpp:44
MoFEM::OpLoopSide
Element used to execute operators on side of the element.
Definition: ForcesAndSourcesCore.hpp:1290
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
MOFEM_INVALID_DATA
@ MOFEM_INVALID_DATA
Definition: definitions.h:36
F
@ F
Definition: free_surface.cpp:394
MoFEM::PetscOptionsGetBool
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
Definition: DeprecatedPetsc.hpp:182
MoFEM::OpPostProcMapInMoab
Post post-proc data at points from hash maps.
Definition: PostProcBrokenMeshInMoabBase.hpp:698