v0.13.0
test_cache_on_entities.cpp
Go to the documentation of this file.
1 /**
2  * \file test_cache_on_entities.cpp
3  * \example test_cache_on_entities.cpp
4  *
5  * Tetsing entities cache. Entities acache is data abstraction enabling user to
6  * pass data between operators, finite elements, and problems, in convinient
7  * way.
8  *
9  * It can be used to store indices, or history variables, and any other data,
10  * which I can not imagine at that moment of time, and how to make use of it
11  * depends on user creativity and imagination.
12  *
13  */
14 
15 /* This file is part of MoFEM.
16  * MoFEM is free software: you can redistribute it and/or modify it under
17  * the terms of the GNU Lesser General Public License as published by the
18  * Free Software Foundation, either version 3 of the License, or (at your
19  * option) any later version.
20  *
21  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
22  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
23  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
24  * License for more details.
25  *
26  * You should have received a copy of the GNU Lesser General Public
27  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>. */
28 
29 #include <MoFEM.hpp>
30 using namespace MoFEM;
31 
32 static char help[] = "...\n\n";
33 
38 
40  MyStorage(VectorInt &data) : globalIndices(data) {}
42 };
43 
44 /**
45  * @brief Operator set cache stored data, in this example, global indices, but
46  * it can be any structure
47  *
48  */
49 struct OpVolumeSet : public VolOp {
50  OpVolumeSet(const std::string &field_name) : VolOp(field_name, OPROW) {}
51  MoFEMErrorCode doWork(int side, EntityType type,
53 
55  MOFEM_LOG_CHANNEL("SYNC");
56  MOFEM_LOG_TAG("SYNC", "OpVolumeSet");
57 
58  // Clear data when start process element
59  if (type == MBVERTEX)
60  entsIndices.clear();
61 
62  // Get field entity pointer
63  auto field_ents = data.getFieldEntities();
64  if (auto e_ptr = field_ents[0]) {
65  // Add indices to global storage
66  entsIndices.push_back(boost::make_shared<MyStorage>(data.getIndices()));
67  // Store pointer to data on entity
68  e_ptr->getWeakStoragePtr() = entsIndices.back();
69 
70  MOFEM_LOG("SYNC", Sev::inform) << "Set " << e_ptr->getEntTypeName() << " "
71  << side << " : " << entsIndices.size();
72 
73  // Check if all works
74  if (auto ptr = e_ptr->getSharedStoragePtr<EntityStorage>()) {
75 
76  if (auto cast_ptr = boost::dynamic_pointer_cast<MyStorage>(ptr))
77  MOFEM_LOG("SYNC", Sev::noisy) << "Cast works";
78  else
79  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Cast not works");
80 
81  } else {
82  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
83  "Pointer to stored data is expected to be set");
84  }
85 
86  if (auto ptr = e_ptr->getSharedStoragePtr<MyStorage>()) {
87  MOFEM_LOG("SYNC", Sev::verbose)
88  << data.getIndices() << " : " << ptr->globalIndices;
89  } else {
90  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
91  "Pointer to stored data is expected to be set");
92  }
93 
94  } else {
95  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
96  "Pointer to entity should be set by finite element");
97  }
98 
100  }
101 
102  using SharedVecInt = boost::shared_ptr<MyStorage>;
103  static std::vector<SharedVecInt>
104  entsIndices; ///< This is global static storage
105 };
106 
107 std::vector<OpVolumeSet::SharedVecInt> OpVolumeSet::entsIndices;
108 
109 /**
110  * @brief Test if cached data can be accessed, and check consistency of data
111  *
112  */
113 struct OpVolumeTest : public VolOp {
114  OpVolumeTest(const std::string &field_name) : VolOp(field_name, OPROW) {}
115  MoFEMErrorCode doWork(int side, EntityType type,
117 
119  MOFEM_LOG_CHANNEL("SYNC");
120  MOFEM_LOG_TAG("SYNC", "OpVolumeTest");
121 
122  // Get pointer to field entities
123  auto field_ents = data.getFieldEntities();
124  if (auto e_ptr = field_ents[0]) {
125 
126  MOFEM_LOG("SYNC", Sev::inform)
127  << "Test " << e_ptr->getEntTypeName() << " " << side;
128 
129  // Check if data are cached on entity, and if code is correct, data should
130  // accessible.
131  if (auto ptr = e_ptr->getSharedStoragePtr<MyStorage>()) {
132 
133  MOFEM_LOG("SYNC", Sev::verbose)
134  << data.getIndices() << " : " << ptr->globalIndices;
135 
136  // Check constancy of data. Stored data are indices, and expected stored
137  // that should be indices, thus difference between should be zero.
138  auto diff = data.getIndices() - ptr->globalIndices;
139  auto dot = inner_prod(diff, diff);
140  if (dot > 0)
141  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID, "Test failed");
142 
143  } else {
144  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
145  "Pointer to stored data is expected to be set");
146  }
147 
148  } else {
149  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
150  "Pointer to entity should be set by finite element");
151  }
152 
154  }
155 };
156 
157 struct OpVolumeAssemble : public VolOp {
158  OpVolumeAssemble(const std::string &field_name, Vec v)
159  : VolOp(field_name, OPROW), V(v) {}
160  MoFEMErrorCode doWork(int side, EntityType type,
163  auto nb_dofs = data.getIndices().size();
164  if (nb_dofs) {
165  std::vector<double> v(nb_dofs, 1);
166  CHKERR VecSetValues<EssentialBcStorage>(V, data, &v[0], ADD_VALUES);
167  }
169  }
170 
171 private:
173 };
174 
175 struct VolRule {
176  int operator()(int, int, int) const { return 1; }
177 };
178 
179 int main(int argc, char *argv[]) {
180 
181  // initialize petsc
182  MoFEM::Core::Initialize(&argc, &argv, (char *)0, help);
183 
184  try {
185 
186  // Create MoAB database
187  moab::Core moab_core;
188  moab::Interface &moab = moab_core;
189 
190  // Create MoFEM database and link it to MoAB
191  MoFEM::Core mofem_core(moab);
192  MoFEM::Interface &m_field = mofem_core;
193 
194  // Register DM Manager
195  DMType dm_name = "DMMOFEM";
196  CHKERR DMRegister_MoFEM(dm_name);
197 
198  // Simple interface
199  Simple *simple_interface;
200  CHKERR m_field.getInterface(simple_interface);
201  {
202  // get options from command line
203  CHKERR simple_interface->getOptions();
204  // load mesh file
205  CHKERR simple_interface->loadFile();
206  // add fields
207  CHKERR simple_interface->addDomainField("FIELD", H1,
209  // set fields order
210  CHKERR simple_interface->setFieldOrder("FIELD", 4);
211  // setup problem
212  CHKERR simple_interface->setUp();
213 
214  // get dm
215  auto dm = simple_interface->getDM();
216  // create elements
217  auto domain_fe = boost::make_shared<VolEle>(m_field);
218  // set integration rule
219  domain_fe->getRuleHook = VolRule();
220 
221  auto get_skin = [&]() {
222  Range ents;
223  CHKERR m_field.get_moab().get_entities_by_dimension(0, 3, ents, true);
224  Skinner skin(&m_field.get_moab());
225  Range skin_ents;
226  CHKERR skin.find_skin(0, ents, false, skin_ents);
227  // filter not owned entities, those are not on boundary
228  ParallelComm *pcomm =
229  ParallelComm::get_pcomm(&m_field.get_moab(), MYPCOMM_INDEX);
230  Range proc_skin;
231  CHKERR pcomm->filter_pstatus(skin_ents,
232  PSTATUS_SHARED | PSTATUS_MULTISHARED,
233  PSTATUS_NOT, -1, &proc_skin);
234  Range adj;
235  for (auto d : {0, 1, 2})
236  CHKERR m_field.get_moab().get_adjacencies(proc_skin, d, false, adj,
237  moab::Interface::UNION);
238  // proc_skin.merge(adj);
239  return proc_skin;
240  };
241 
242  auto get_mark_skin_dofs = [&](Range &&skin) {
243  auto problem_manager = m_field.getInterface<ProblemsManager>();
244  auto marker_ptr = boost::make_shared<std::vector<unsigned char>>();
245  problem_manager->markDofs(simple_interface->getProblemName(), ROW, skin,
246  *marker_ptr);
247  return marker_ptr;
248  };
249 
252  auto mark_dofs = get_mark_skin_dofs(get_skin());
253 
254  // set operator to the volume elements
255  domain_fe->getOpPtrVector().push_back(new OpVolumeSet("FIELD"));
256  domain_fe->getOpPtrVector().push_back(new OpVolumeTest("FIELD"));
257 
258  domain_fe->getOpPtrVector().push_back(
259  new OpSetBc("FIELD", true, mark_dofs));
260  domain_fe->getOpPtrVector().push_back(new OpVolumeAssemble("FIELD", v));
261  domain_fe->getOpPtrVector().push_back(new OpUnSetBc("FIELD"));
262 
263  // make integration in volume (here real calculations starts)
264  CHKERR DMoFEMLoopFiniteElements(dm, simple_interface->getDomainFEName(),
265  domain_fe);
266  OpVolumeSet::entsIndices.clear();
267 
268  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
269  CHKERR VecAssemblyBegin(v);
270  CHKERR VecAssemblyEnd(v);
271  CHKERR VecGhostUpdateBegin(v, ADD_VALUES, SCATTER_REVERSE);
272  CHKERR VecGhostUpdateEnd(v, ADD_VALUES, SCATTER_REVERSE);
273  CHKERR VecGhostUpdateBegin(v, INSERT_VALUES, SCATTER_FORWARD);
274  CHKERR VecGhostUpdateEnd(v, INSERT_VALUES, SCATTER_FORWARD);
275 
276  CHKERR DMoFEMMeshToLocalVector(simple_interface->getDM(), v,
277  INSERT_VALUES, SCATTER_REVERSE);
278 
279  // CHKERR VecView(v,PETSC_VIEWER_STDOUT_WORLD);
280 
281  const MoFEM::Problem *problem_ptr;
282  CHKERR DMMoFEMGetProblemPtr(dm, &problem_ptr);
283 
284  double *array;
285  CHKERR VecGetArray(v, &array);
286  for (auto dof : *(problem_ptr->getNumeredRowDofsPtr())) {
287  auto loc_idx = dof->getPetscLocalDofIdx();
288 
289  if (dof->getFieldData() != array[loc_idx])
290  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
291  "Data inconsistency");
292 
293  if ((*mark_dofs)[loc_idx]) {
294  if (array[loc_idx] != 0)
295  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
296  "Dof assembled, but it should be not");
297  } else {
298  if (array[loc_idx] == 0)
299  SETERRQ(PETSC_COMM_SELF, MOFEM_ATOM_TEST_INVALID,
300  "Dof not assembled, but it should be");
301  }
302  }
303 
304  CHKERR VecRestoreArray(v, &array);
305  }
306  }
307  CATCH_ERRORS;
308 
309  // finish work cleaning memory, getting statistics, etc.
311 
312  return 0;
313 }
ForcesAndSourcesCore::UserDataOperator UserDataOperator
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:348
MoFEM::FaceElementForcesAndSourcesCore FaceEle
@ ROW
Definition: definitions.h:136
#define CATCH_ERRORS
Catch errors.
Definition: definitions.h:385
@ AINSWORTH_LEGENDRE_BASE
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:73
@ H1
continuous field
Definition: definitions.h:98
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:228
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:359
@ MOFEM_ATOM_TEST_INVALID
Definition: definitions.h:53
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:44
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:429
#define CHKERR
Inline error check.
Definition: definitions.h:548
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:481
PetscErrorCode DMMoFEMGetProblemPtr(DM dm, const MoFEM::Problem **problem_ptr)
Get pointer to problem data structure.
Definition: DMMMoFEM.cpp:384
PetscErrorCode DMRegister_MoFEM(const char sname[])
Register MoFEM problem.
Definition: DMMMoFEM.cpp:59
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition: DMMMoFEM.cpp:544
PetscErrorCode DMCreateGlobalVector_MoFEM(DM dm, Vec *g)
DMShellSetCreateGlobalVector.
Definition: DMMMoFEM.cpp:1105
FaceElementForcesAndSourcesCoreSwitch< 0 > FaceElementForcesAndSourcesCore
Face finite element default.
VolumeElementForcesAndSourcesCoreSwitch< 0 > VolumeElementForcesAndSourcesCore
Volume finite element default.
MoFEMErrorCode VecSetValues< EssentialBcStorage >(Vec V, const EntitiesFieldData::EntData &data, const double *ptr, InsertMode iora)
Set values to vector in operator.
#define MOFEM_LOG(channel, severity)
Log.
Definition: LogManager.hpp:311
#define MOFEM_LOG_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:342
#define MOFEM_LOG_CHANNEL(channel)
Set and reset channel.
Definition: LogManager.hpp:287
double v
phase velocity of light in medium (cm/ns)
const FTensor::Tensor2< T, Dim, Dim > Vec
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T * > &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:67
UBlasVector< int > VectorInt
Definition: Types.hpp:78
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
CoreTmp< 0 > Core
Definition: Core.hpp:1096
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1965
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition: Core.hpp:92
static MoFEMErrorCode Initialize(int *argc, char ***args, const char file[], const char help[])
Initializes the MoFEM database PETSc, MOAB and MPI.
Definition: Core.cpp:85
static MoFEMErrorCode Finalize()
Checks for options to be called at the conclusion of the program.
Definition: Core.cpp:125
Deprecated interface functions.
Data on single entity (This is passed as argument to DataOperator::doWork)
const VectorFieldEntities & getFieldEntities() const
get field entities
const VectorInt & getIndices() const
Get global indices of dofs on entity.
Set indices on entities on finite element.
keeps basic data about problem
auto & getNumeredRowDofsPtr() const
get access to numeredRowDofsPtr storing DOFs on rows
Problem manager is used to build and partition problems.
MoFEMErrorCode markDofs(const std::string problem_name, RowColData rc, const enum MarkOP op, const Range ents, std::vector< unsigned char > &marker) const
Create vector with marked indices.
Simple interface for fast problem set-up.
Definition: Simple.hpp:33
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:293
MoFEMErrorCode getOptions()
get options
Definition: Simple.cpp:216
MoFEMErrorCode getDM(DM *dm)
Get DM.
Definition: Simple.cpp:715
MoFEMErrorCode loadFile(const std::string options, const std::string mesh_file_name)
Load mesh file.
Definition: Simple.cpp:230
MoFEMErrorCode setFieldOrder(const std::string field_name, const int order, const Range *ents=NULL)
Set field order.
Definition: Simple.cpp:508
MoFEMErrorCode setUp(const PetscBool is_partitioned=PETSC_TRUE)
Setup problem.
Definition: Simple.cpp:668
const std::string getProblemName() const
Get the Problem Name.
Definition: Simple.hpp:340
const std::string getDomainFEName() const
Get the Domain FE Name.
Definition: Simple.hpp:319
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
VolumeElementForcesAndSourcesCoreBase::UserDataOperator UserDataOperator
MyStorage(VectorInt &data)
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpVolumeAssemble(const std::string &field_name, Vec v)
Operator set cache stored data, in this example, global indices, but it can be any structure.
boost::shared_ptr< MyStorage > SharedVecInt
OpVolumeSet(const std::string &field_name)
static std::vector< SharedVecInt > entsIndices
This is global static storage.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
Test if cached data can be accessed, and check consistency of data.
MoFEMErrorCode doWork(int side, EntityType type, EntitiesFieldData::EntData &data)
Operator for linear form, usually to calculate values on right hand side.
OpVolumeTest(const std::string &field_name)
Set integration rule.
int operator()(int, int, int) const
int main(int argc, char *argv[])
static char help[]
FaceEle::UserDataOperator FaceOp