v0.11.1
FieldEvaluator.cpp
Go to the documentation of this file.
1 /** \file FieldEvaluator.cpp
2  * \brief TetGen interface for meshing/re-meshing and on the fly mesh creation
3  *
4  */
5 
6 /**
7  * MoFEM is free software: you can redistribute it and/or modify it under
8  * the terms of the GNU Lesser General Public License as published by the
9  * Free Software Foundation, either version 3 of the License, or (at your
10  * option) any later version.
11  *
12  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
13  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  * License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
19  */
20 
21 namespace MoFEM {
22 
24  : cOre(const_cast<MoFEM::Core &>(core)) {
25 
26  if (!LogManager::checkIfChannelExist("FieldEvaluatorWorld")) {
27  auto core_log = logging::core::get();
28 
30  "FieldEvaluatorWorld"));
32  "FieldEvaluatorSync"));
34  "FieldEvaluatorSelf"));
35 
36  LogManager::setLog("FieldEvaluatorWorld");
37  LogManager::setLog("FieldEvaluatorSync");
38  LogManager::setLog("FieldEvaluatorSelf");
39 
40  MOFEM_LOG_TAG("FieldEvaluatorWorld", "FieldEvaluator");
41  MOFEM_LOG_TAG("FieldEvaluatorSync", "FieldEvaluator");
42  MOFEM_LOG_TAG("FieldEvaluatorSelf", "FieldEvaluator");
43  }
44 
45  MOFEM_LOG("FieldEvaluatorWorld", Sev::noisy) << "Field evaluator intreface";
46 }
47 
50  UnknownInterface **iface) const {
52  *iface = NULL;
53  if (uuid == IDD_MOFEMFieldEvaluator) {
54  *iface = const_cast<FieldEvaluatorInterface *>(this);
56  }
57  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
59 }
60 
62 FieldEvaluatorInterface::buildTree3D(boost::shared_ptr<SetPtsData> spd_ptr,
63  const std::string finite_element) {
64  CoreInterface &m_field = cOre;
66  EntityHandle fe_meshset;
67  fe_meshset = m_field.get_finite_element_meshset(finite_element);
68  Range entities_3d;
69  CHKERR m_field.get_moab().get_entities_by_dimension(fe_meshset, 3,
70  entities_3d, true);
71  spd_ptr->treePtr.reset(new AdaptiveKDTree(&m_field.get_moab()));
72  CHKERR spd_ptr->treePtr->build_tree(entities_3d, &spd_ptr->rooTreeSet);
74 }
75 
77 operator()(int order_row, int order_col, int order_data) {
79 
80  if (auto data_ptr = dataPtr.lock()) {
81 
82  decltype(data_ptr->nbEvalPoints) nbEvalPoints = data_ptr->nbEvalPoints;
83  decltype(data_ptr->verb) verb = data_ptr->verb;
84 
85  decltype(data_ptr->shapeFunctions) &shapeFunctions =
86  data_ptr->shapeFunctions;
87  decltype(data_ptr->evalPointEntityHandle) &evalPointEntityHandle =
88  data_ptr->evalPointEntityHandle;
89 
90  if (auto fe_ptr = data_ptr->feMethodPtr.lock()) {
91 
93  static_cast<VolumeElementForcesAndSourcesCore &>(*fe_ptr);
94 
95  MatrixDouble &gauss_pts = fe.gaussPts;
96  gauss_pts.resize(4, nbEvalPoints, false);
97  gauss_pts.clear();
98 
100  &shapeFunctions(0, 0), &shapeFunctions(0, 1), &shapeFunctions(0, 2),
101  &shapeFunctions(0, 3)};
103  &gauss_pts(0, 0), &gauss_pts(1, 0), &gauss_pts(2, 0)};
104 
105  int nb_gauss_pts = 0;
106  for (int nn = 0; nn != nbEvalPoints; ++nn) {
107  if (evalPointEntityHandle[nn] ==
108  fe.numeredEntFiniteElementPtr->getEnt()) {
109  for (const int i : {0, 1, 2}) {
110  t_gauss_pts(i) = t_shape(i + 1);
111  gauss_pts(3, nb_gauss_pts) = nn;
112  }
113  ++t_gauss_pts;
114  ++nb_gauss_pts;
115  }
116  ++t_shape;
117  }
118 
119  if (verb >= VERY_NOISY)
120  std::cout << "nbEvalOPoints / nbGaussPts: " << nbEvalPoints << " / "
121  << nb_gauss_pts << std::endl;
122  gauss_pts.resize(4, nb_gauss_pts, true);
123 
124  if (verb >= VERY_NOISY)
125  std::cout << "gauss pts: " << gauss_pts << std::endl;
126 
127  } else
128  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
129  "Pointer to element does not exists");
130 
131  } else
132  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
133  "Pointer to data does not exists");
134 
136 }
137 
139  const double *const point, const double distance, const std::string problem,
140  const std::string finite_element, boost::shared_ptr<SetPtsData> data_ptr,
141  int lower_rank, int upper_rank, boost::shared_ptr<CacheTuple> cache_ptr,
142  MoFEMTypes bh, VERBOSITY_LEVELS verb) {
143  CoreInterface &m_field = cOre;
145 
146  std::vector<EntityHandle> leafs_out;
147  CHKERR data_ptr->treePtr->distance_search(point, distance, leafs_out);
148  Range tree_ents;
149  for (auto lit : leafs_out)
150  CHKERR m_field.get_moab().get_entities_by_dimension(lit, 3, tree_ents,
151  true);
152 
153  if (verb >= VERY_NOISY)
154  std::cout << "tree entities: " << tree_ents << endl;
155 
156  data_ptr->evalPointEntityHandle.resize(data_ptr->nbEvalPoints);
157  std::fill(data_ptr->evalPointEntityHandle.begin(),
158  data_ptr->evalPointEntityHandle.end(), 0);
159 
160  for (auto tet : tree_ents) {
161 
162  const EntityHandle *conn;
163  int num_nodes;
164  CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes, true);
165  std::array<double, 12> coords;
166  CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
167 
168  for (int n = 0; n != data_ptr->nbEvalPoints; ++n) {
169 
170  std::array<double, 3> local_coords;
172  coords.data(), &data_ptr->evalPoints[3 * n], 1, local_coords.data());
173 
174  std::array<double, 4> shape;
175  CHKERR Tools::shapeFunMBTET<3>(shape.data(), &local_coords[0],
176  &local_coords[1], &local_coords[2], 1);
177 
178  const double eps = data_ptr->eps;
179  if (shape[0] >= 0 - eps && shape[0] <= 1 + eps &&
180 
181  shape[1] >= 0 - eps && shape[1] <= 1 + eps &&
182 
183  shape[2] >= 0 - eps && shape[2] <= 1 + eps &&
184 
185  shape[3] >= 0 - eps && shape[3] <= 1 + eps) {
186 
187  std::copy(shape.begin(), shape.end(), &data_ptr->shapeFunctions(n, 0));
188  std::copy(local_coords.begin(), local_coords.end(),
189  &data_ptr->localCoords(n, 0));
190  data_ptr->evalPointEntityHandle[n] = tet;
191  }
192  }
193  }
194 
195  const Problem *prb_ptr;
196  CHKERR m_field.get_problem(problem, &prb_ptr);
197  boost::shared_ptr<NumeredEntFiniteElement_multiIndex> numered_fes(
199 
200  Range in_tets;
201  in_tets.insert_list(data_ptr->evalPointEntityHandle.begin(),
202  data_ptr->evalPointEntityHandle.end());
203  in_tets = in_tets.subset_by_dimension(3);
204 
205  if (verb >= VERY_NOISY)
206  std::cout << "in tets: " << in_tets << endl;
207 
208  for (auto peit = in_tets.pair_begin(); peit != in_tets.pair_end(); ++peit) {
209  auto lo =
211  .lower_bound(boost::make_tuple(finite_element, peit->first));
212  auto hi =
214  .upper_bound(boost::make_tuple(finite_element, peit->second));
215  numered_fes->insert(lo, hi);
216 
217  if (verb >= VERY_NOISY)
218  std::cout << "numered elements:" << std::endl;
219  for (; lo != hi; ++lo)
220  if (verb >= VERY_NOISY)
221  std::cout << **lo << endl;
222  }
223  if (verb >= VERY_NOISY)
224  std::cout << std::endl;
225 
226  if (auto fe_ptr = data_ptr->feMethodPtr.lock()) {
227 
228  MOFEM_LOG_C("FieldEvaluatorSync", Sev::verbose,
229  "Number elements %d to evaluate at proc %d",
230  numered_fes->size(), m_field.get_comm_rank());
231  MOFEM_LOG_SYNCHRONISE(m_field.get_comm());
232 
233  if(!cache_ptr) {
234  cache_ptr = boost::make_shared<CacheTuple>();
235  CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
236 
237  MOFEM_LOG("FieldEvaluatorSync", Sev::noisy)
238  << "If you call that function many times in the loop consider to set "
239  "cache_ptr outside of the loop. Otherwise code can be slow.";
240  }
241 
242  CHKERR m_field.loop_finite_elements(prb_ptr, finite_element, *fe_ptr,
243  lower_rank, upper_rank, numered_fes, bh,
244  cache_ptr, verb);
245 
246  } else
247  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
248  "Pointer to element does not exists");
249 
251 }
252 
253 } // namespace MoFEM
static Index< 'n', 3 > n
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
Definition: LogManager.hpp:349
#define MOFEM_LOG_C(channel, severity, format,...)
Definition: LogManager.hpp:315
VERBOSITY_LEVELS
Verbosity levels.
Definition: definitions.h:275
@ VERY_NOISY
Definition: definitions.h:281
MoFEMTypes
Those types control how functions respond on arguments, f.e. error handling.
Definition: definitions.h:189
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:516
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:415
@ MOFEM_DATA_INCONSISTENCY
Definition: definitions.h:123
#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
multi_index_container< boost::shared_ptr< NumeredEntFiniteElement >, indexed_by< ordered_unique< tag< Unique_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_EntFiniteElement, UId, &NumeredEntFiniteElement::getLocalUniqueId > >, ordered_non_unique< tag< Part_mi_tag >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > >, ordered_non_unique< tag< FiniteElement_name_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef > >, ordered_non_unique< tag< Ent_mi_tag >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > >, ordered_non_unique< tag< Composite_Name_And_Ent_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, const_mem_fun< NumeredEntFiniteElement::interface_type_RefEntity, EntityHandle, &NumeredEntFiniteElement::getEnt > > >, ordered_non_unique< tag< Composite_Name_And_Part_mi_tag >, composite_key< NumeredEntFiniteElement, const_mem_fun< NumeredEntFiniteElement::interface_type_FiniteElement, boost::string_ref, &NumeredEntFiniteElement::getNameRef >, member< NumeredEntFiniteElement, unsigned int, &NumeredEntFiniteElement::part > > > > > NumeredEntFiniteElement_multiIndex
MultiIndex for entities for NumeredEntFiniteElement.
virtual const Problem * get_problem(const std::string problem_name) const =0
Get the problem object.
virtual EntityHandle get_finite_element_meshset(const std::string &name) const =0
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_TAG(channel, tag)
Tag channel.
Definition: LogManager.hpp:343
virtual MoFEMErrorCode loop_finite_elements(const std::string &problem_name, const std::string &fe_name, FEMethod &method, boost::shared_ptr< NumeredEntFiniteElement_multiIndex > fe_ptr=nullptr, MoFEMTypes bh=MF_EXIST, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr(), int verb=DEFAULT_VERBOSITY)=0
Make a loop over finite elements.
FTensor::Index< 'i', 3 > i
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
static const MOFEMuuid IDD_MOFEMFieldEvaluator
virtual MoFEMErrorCode cache_problem_entities(const std::string prb_name, CacheTupleWeakPtr cache_ptr)=0
Cache variables.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
virtual int get_comm_rank() const =0
Core (interface) class.
Definition: Core.hpp:77
boost::shared_ptr< const NumeredEntFiniteElement > numeredEntFiniteElementPtr
MoFEMErrorCode operator()(int order_row, int order_col, int order_data)
boost::weak_ptr< SetPtsData > dataPtr
Field evaluator interface.
MoFEMErrorCode buildTree3D(boost::shared_ptr< SetPtsData > spd_ptr, const std::string finite_element)
Build spatial tree.
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
FieldEvaluatorInterface(const MoFEM::Core &core)
MoFEMErrorCode evalFEAtThePoint3D(const double *const point, const double distance, const std::string problem, const std::string finite_element, boost::shared_ptr< SetPtsData > data_ptr, int lower_rank, int upper_rank, boost::shared_ptr< CacheTuple > cache_ptr, MoFEMTypes bh=MF_EXIST, VERBOSITY_LEVELS verb=QUIET)
Evaluate field at artbitray position.
MatrixDouble gaussPts
Matrix of integration points.
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 > getStrmWorld()
Get the strm world object.
Definition: LogManager.cpp:327
static boost::shared_ptr< std::ostream > getStrmSync()
Get the strm sync object.
Definition: LogManager.cpp:331
static bool checkIfChannelExist(const std::string channel)
Check if channel exist.
Definition: LogManager.cpp:378
static boost::shared_ptr< std::ostream > getStrmSelf()
Get the strm self object.
Definition: LogManager.cpp:323
MoFEM interface unique ID.
keeps basic data about problem
std::string getName() const
boost::shared_ptr< NumeredEntFiniteElement_multiIndex > numeredFiniteElementsPtr
store finite elements
static MoFEMErrorCode getLocalCoordinatesOnReferenceFourNodeTet(const double *elem_coords, const double *glob_coords, const int nb_nodes, double *local_coords)
Get the Local Coordinates On Reference Four Node Tet object.
Definition: Tools.cpp:97
base class for all interface classes