v0.11.1
continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
Go to the documentation of this file.
1 /**
2  * \file continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
3  * \ingroup mofem_simple_interface
4  \example continuity_check_on_skeleton_with_simple_2d_for_h1.cpp
5  *
6  * \brief Integration on skeleton for 2d
7  *
8  * Teting integration on skeleton and checking of continuity of hcurl space on
9  * edges.
10  */
11 
12 /* This file is part of MoFEM.
13  * MoFEM is free software: you can redistribute it and/or modify it under
14  * the terms of the GNU Lesser General Public License as published by the
15  * Free Software Foundation, either version 3 of the License, or (at your
16  * option) any later version.
17  *
18  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
19  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
21  * License for more details.
22  *
23  * You should have received a copy of the GNU Lesser General Public
24  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
25 
26 #include <MoFEM.hpp>
27 using namespace MoFEM;
28 
29 static char help[] = "...\n\n";
30 
32  FaceElementForcesAndSourcesCore::NO_HO_GEOMETRY>;
33 
36 
39 
40 struct CommonData {
44 };
45 
46 struct SkeletonFE : public EdgeEleOp {
47 
48  struct OpFaceSide : public FaceEleOnSideOp {
49 
51  OpFaceSide(CommonData &elem_data)
52  : FaceEleOnSideOp("FIELD", UserDataOperator::OPROW),
53  elemData(elem_data) {}
54 
55  MoFEMErrorCode doWork(int side, EntityType type,
58 
59  if (type == MBEDGE && side == getEdgeSideNumber()) {
60 
61  MatrixDouble diff = getCoordsAtGaussPts() - getEdgeCoordsAtGaussPts();
62 
63  const double eps = 1e-12;
64  if (norm_inf(diff) > eps) {
65  MOFEM_LOG_ATTRIBUTES("ATOM",
68  MOFEM_LOG("ATOM", Sev::error)
69  << "Quad coords: " << getCoordsAtGaussPts();
70  MOFEM_LOG("ATOM", Sev::error)
71  << "Edge coords: " << getEdgeCoordsAtGaussPts();
72  MOFEM_LOG("ATOM", Sev::error) << "Difference: " << diff;
73  SETERRQ(PETSC_COMM_WORLD, MOFEM_ATOM_TEST_INVALID,
74  "Coordinates at integration pts are different");
75  }
76 
77  const size_t nb_dofs = data.getN().size2();
78  const size_t nb_integration_pts = data.getN().size1();
79 
80  auto t_base = data.getFTensor0N();
81  MatrixDouble *ptr_dot_elem_data = nullptr;
82  if (getFEMethod()->nInTheLoop == 0)
83  ptr_dot_elem_data = &elemData.dotEleLeft;
84  else
85  ptr_dot_elem_data = &elemData.dotEleRight;
86  MatrixDouble &dot_elem_data = *ptr_dot_elem_data;
87  dot_elem_data.resize(nb_integration_pts, nb_dofs, false);
88 
89  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
90  for (size_t bb = 0; bb != nb_dofs; ++bb) {
91  dot_elem_data(gg, bb) = t_base;
92  ++t_base;
93  }
94  }
95  }
97  }
98  };
99 
102 
103  SkeletonFE(MoFEM::Interface &m_field, CommonData &elem_data)
104  : EdgeEle::UserDataOperator("FIELD", UserDataOperator::OPROW),
105  faceSideFe(m_field), elemData(elem_data) {
106  faceSideFe.getOpPtrVector().push_back(new SkeletonFE::OpFaceSide(elemData));
107  }
108 
109  MoFEMErrorCode doWork(int side, EntityType type,
111 
113  if (type == MBEDGE) {
114 
115  const size_t nb_dofs = data.getN().size2();
116  const size_t nb_integration_pts = data.getN().size1();
117 
118  MOFEM_LOG("ATOM", Sev::noisy)
119  << "Cords at integration points" << getCoordsAtGaussPts();
120 
121  auto t_base = data.getFTensor0N();
122  elemData.dotEdge.resize(nb_integration_pts, nb_dofs, false);
123  elemData.dotEleLeft.resize(nb_integration_pts, 0, false);
124  elemData.dotEleRight.resize(nb_integration_pts, 0, false);
125 
126  for (size_t gg = 0; gg != nb_integration_pts; ++gg) {
127  for (size_t bb = 0; bb != nb_dofs; ++bb) {
128  elemData.dotEdge(gg, bb) = t_base;
129  ++t_base;
130  }
131  }
132 
133  CHKERR loopSideFaces("dFE", faceSideFe);
134 
135  auto check_continuity_of_base = [&](auto &vol_dot_data) {
137 
138  if (vol_dot_data.size1() != elemData.dotEdge.size1())
139  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
140  "Inconsistent number of integration points");
141 
142  if (vol_dot_data.size2() != elemData.dotEdge.size2())
143  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
144  "Inconsistent number of base functions");
145  const double eps = 1e-12;
146  for (size_t gg = 0; gg != vol_dot_data.size1(); ++gg)
147  for (size_t bb = 0; bb != vol_dot_data.size2(); ++bb) {
148  const double error =
149  std::abs(vol_dot_data(gg, bb) - elemData.dotEdge(gg, bb));
150  if (error > eps)
151  SETERRQ4(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
152  "Inconsistency (%d, %d) %3.4e != %3.4e", gg, bb,
153  vol_dot_data(gg, bb), elemData.dotEdge(gg, bb));
154  else
155  MOFEM_LOG("ATOM", Sev::noisy) << "Ok";
156  }
158  };
159 
160  if (elemData.dotEleLeft.size2() != 0)
161  CHKERR check_continuity_of_base(elemData.dotEleLeft);
162  else if (elemData.dotEleRight.size2() != 0)
163  CHKERR check_continuity_of_base(elemData.dotEleRight);
164  }
166  }
167 };
168 
169 int main(int argc, char *argv[]) {
170 
171  // initialize petsc
172  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
173 
174  try {
175 
176  // Create MoAB database
177  moab::Core moab_core;
178  moab::Interface &moab = moab_core;
179 
180  // Create MoFEM database and link it to MoAB
181  MoFEM::Core mofem_core(moab);
182  MoFEM::Interface &m_field = mofem_core;
183 
184  auto core_log = logging::core::get();
185  core_log->add_sink(
187  LogManager::setLog("ATOM");
188 
189  // Register DM Manager
190  DMType dm_name = "DMMOFEM";
191  CHKERR DMRegister_MoFEM(dm_name);
192 
193  // Simple interface
194  Simple *simple_interface;
195  CHKERR m_field.getInterface(simple_interface);
196  {
197  // get options from command line
198  CHKERR simple_interface->getOptions();
199  // load mesh file
200  CHKERR simple_interface->loadFile("");
201 
202  auto get_base = []() -> FieldApproximationBase {
203  enum bases { AINSWORTH, DEMKOWICZ, LASTBASEOP };
204  const char *list_bases[] = {"ainsworth", "demkowicz"};
205  PetscBool flg;
206  PetscInt choice_base_value = AINSWORTH;
207  CHKERR PetscOptionsGetEList(PETSC_NULL, NULL, "-base", list_bases,
208  LASTBASEOP, &choice_base_value, &flg);
209  if (flg == PETSC_TRUE) {
211  if (choice_base_value == AINSWORTH)
213  else if (choice_base_value == DEMKOWICZ)
214  base = DEMKOWICZ_JACOBI_BASE;
215  return base;
216  }
217  return LASTBASE;
218  };
219 
220  // add fields
221  auto base = get_base();
222  CHKERR simple_interface->addDomainField("FIELD", H1, base, 1);
223  // CHKERR simple_interface->addDomainField("TEST_FIELD", L2,
224  // AINSWORTH_LEGENDRE_BASE, 1);
225  CHKERR simple_interface->addSkeletonField("FIELD", H1, base, 1);
226  // set fields order
227  CHKERR simple_interface->setFieldOrder("FIELD", 2);
228  // CHKERR simple_interface->setFieldOrder("TEST_FIELD", 1);
229  // setup problem
230  CHKERR simple_interface->setUp();
231  // get dm
232  auto dm = simple_interface->getDM();
233 
234  // create elements
235  CommonData elem_data;
236  boost::shared_ptr<EdgeEle> skeleton_fe =
237  boost::shared_ptr<EdgeEle>(new EdgeEle(m_field));
238  skeleton_fe->getOpPtrVector().push_back(
239  new SkeletonFE(m_field, elem_data));
240 
241  // iterate skeleton finite elements
242  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getSkeletonFEName(),
243  skeleton_fe);
244  }
245  }
246  CATCH_ERRORS;
247 
248  // finish work cleaning memory, getting statistics, etc.
250 
251  return 0;
252 }
int main(int argc, char *argv[])
MoFEM::EdgeElementForcesAndSourcesCoreSwitch< EdgeElementForcesAndSourcesCore::NO_HO_GEOMETRY > EdgeEle
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:441
FieldApproximationBase
approximation base
Definition: definitions.h:150
@ LASTBASE
Definition: definitions.h:161
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
@ DEMKOWICZ_JACOBI_BASE
Definition: definitions.h:158
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
@ H1
continuous field
Definition: definitions.h:177
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:132
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:485
#define CHKERR
Inline error check.
Definition: definitions.h:604
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:509
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method)
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:507
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:48
static LoggerType & setLog(const std::string channel)
Set ans resset chanel logger.
Definition: LogManager.cpp:368
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:312
#define MOFEM_LOG_ATTRIBUTES(channel, bit)
Add attributes to channel.
Definition: LogManager.hpp:300
#define MOFEM_LOG_FUNCTION()
Set scope.
Definition: LogManager.hpp:329
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:76
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1140
PetscErrorCode PetscOptionsGetEList(PetscOptions *, const char pre[], const char name[], const char *const *list, PetscInt next, PetscInt *value, PetscBool *set)
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1953
Core (interface) class.
Definition: Core.hpp:77
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:60
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:100
Data on single entity (This is passed as argument to DataOperator::doWork)
MatrixDouble & getN(const FieldApproximationBase base)
get base functions this return matrix (nb. of rows is equal to nb. of Gauss pts, nb....
FTensor::Tensor0< FTensor::PackPtr< double *, 1 > > getFTensor0N(const FieldApproximationBase base)
Get base function as Tensor0.
Deprecated interface functions.
Data operator to do calculations at integration points.
boost::ptr_vector< UserDataOperator > & getOpPtrVector()
Use to push back operator for row operator.
static boost::shared_ptr< SinkType > createSink(boost::shared_ptr< std::ostream > stream_ptr, std::string comm_filter)
Create a sink object.
Definition: LogManager.cpp:283
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:323
Simple interface for fast problem set-up.
Definition: Simple.hpp:36
const std::string getSkeletonFEName() const
Get the Skeleton FE Name.
Definition: Simple.hpp:322
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:269
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:198
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:713
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:212
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:482
MoFEMErrorCode addSkeletonField(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 skeleton.
Definition: Simple.cpp:305
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:679
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
SkeletonFE(MoFEM::Interface &m_field, CommonData &elem_data)
MoFEMErrorCode doWork(int side, EntityType type, DataForcesAndSourcesCore::EntData &data)