Evaluate field at artbitray position.
std::array<double, 3> point = {0, 0, 0};
const double dist = 0.3;
std::array<double, 6> eval_points = {-1., -1., -1., 1., 1., 1. };
using VolEle = VolumeElementForcesAndSourcesCore;
auto data_ptr =
getData<VolEle>(point.data(), point.size/3);
if(auto vol_ele = data_ptr->feMethod.lock()) {
vol_ele->getOpPtrVector().push_back(
new MyOp());
auto cache_ptr = boost::make_shared<CacheTuple>();
CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
CHKERR m_field.getInterface<FieldEvaluatorInterface>()
"FINITE_ELEMENT_NAME",
data_ptr, m_field.get_comm_rank(),
m_field.get_comm_rank(), cache_ptr);
}
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.
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.
Volume finite element base.
Operator used to check consistency between local coordinates and global cooridnates for integrated po...
165 {
166 CoreInterface &m_field =
cOre;
168
169 std::vector<EntityHandle> leafs_out;
170 CHKERR data_ptr->treePtr->distance_search(point, distance, leafs_out);
172 for (auto lit : leafs_out)
173 CHKERR m_field.get_moab().get_entities_by_dimension(lit,
D, tree_ents,
174 true);
176 MOFEM_LOG(
"SELF", Sev::noisy) <<
"tree entities: " << tree_ents;
177
178 data_ptr->evalPointEntityHandle.resize(data_ptr->nbEvalPoints);
179 std::fill(data_ptr->evalPointEntityHandle.begin(),
180 data_ptr->evalPointEntityHandle.end(), 0);
181
182 std::vector<double> local_coords;
183 std::vector<double> shape;
184
185 auto nb_eval_points = data_ptr->nbEvalPoints;
186
187 for (auto tet : tree_ents) {
189 int num_nodes;
190 CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes,
true);
191
192 if constexpr (
D == 3) {
193
194 local_coords.resize(3 * nb_eval_points);
195 shape.resize(4 * nb_eval_points);
196
197 std::array<double, 12> coords;
198 CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
199
201 coords.data(), &data_ptr->evalPoints[0], nb_eval_points,
202 &local_coords[0]);
203 CHKERR Tools::shapeFunMBTET<3>(&shape[0], &local_coords[0],
204 &local_coords[1], &local_coords[2],
205 nb_eval_points);
206
209 &shape[0], &shape[1], &shape[2], &shape[3]};
211 &data_ptr->shapeFunctions(0, 0), &data_ptr->shapeFunctions(0, 1),
212 &data_ptr->shapeFunctions(0, 2), &data_ptr->shapeFunctions(0, 3)};
213
216 &local_coords[0], &local_coords[1], &local_coords[2]};
218 &(data_ptr->localCoords(0, 0)), &(data_ptr->localCoords(0, 1)),
219 &(data_ptr->localCoords(0, 2))};
220
221 for (
int n = 0;
n != nb_eval_points; ++
n) {
222
223 const double eps = data_ptr->eps;
224 if (t_shape(0) >= 0 -
eps && t_shape(0) <= 1 +
eps &&
225
226 t_shape(1) >= 0 -
eps && t_shape(1) <= 1 +
eps &&
227
228 t_shape(2) >= 0 -
eps && t_shape(2) <= 1 +
eps &&
229
230 t_shape(3) >= 0 -
eps && t_shape(3) <= 1 +
eps) {
231
232 data_ptr->evalPointEntityHandle[
n] = tet;
233 t_shape_data(i4) = t_shape(i4);
234 t_local_data(j3) = t_local(j3);
235 }
236
237 ++t_shape;
238 ++t_shape_data;
239 ++t_local;
240 ++t_local_data;
241 }
242 }
243
244 if constexpr (
D == 2) {
245
246 local_coords.resize(2 * nb_eval_points);
247 shape.resize(3 * nb_eval_points);
248
249 std::array<double, 9> coords;
250 CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
251
253 coords.data(), &data_ptr->evalPoints[0], nb_eval_points,
254 &local_coords[0]);
255 CHKERR Tools::shapeFunMBTRI<2>(&shape[0], &local_coords[0],
256 &local_coords[1], nb_eval_points);
257
260 &shape[0], &shape[1], &shape[2]};
262 &data_ptr->shapeFunctions(0, 0), &data_ptr->shapeFunctions(0, 1),
263 &data_ptr->shapeFunctions(0, 2)};
264
267 &local_coords[0], &local_coords[1]};
268 data_ptr->localCoords.resize(nb_eval_points, 2);
270 &(data_ptr->localCoords(0, 0)), &(data_ptr->localCoords(0, 1))};
271
272 for (
int n = 0;
n != nb_eval_points; ++
n) {
273
274 const double eps = data_ptr->eps;
275 if (t_shape(0) >= 0 -
eps && t_shape(0) <= 1 +
eps &&
276
277 t_shape(1) >= 0 -
eps && t_shape(1) <= 1 +
eps &&
278
279 t_shape(2) >= 0 -
eps && t_shape(2) <= 1 +
eps) {
280
281 data_ptr->evalPointEntityHandle[
n] = tet;
282 t_shape_data(i3) = t_shape(i3);
283 t_local_data(j2) = t_local(j2);
284 }
285
286 ++t_shape;
287 ++t_shape_data;
288 ++t_local;
289 ++t_local_data;
290 }
291 }
292 }
293
294 const Problem *prb_ptr;
295 CHKERR m_field.get_problem(problem, &prb_ptr);
296 boost::shared_ptr<NumeredEntFiniteElement_multiIndex> numered_fes(
298
300 in_tets.insert_list(data_ptr->evalPointEntityHandle.begin(),
301 data_ptr->evalPointEntityHandle.end());
302 in_tets = in_tets.subset_by_dimension(
D);
303
305 MOFEM_LOG(
"SELF", Sev::noisy) <<
"in tets: " << in_tets << endl;
306
307 auto fe_miit =
308 m_field.get_finite_elements()->get<FiniteElement_name_mi_tag>().find(
309 finite_element);
310 if (fe_miit !=
311 m_field.get_finite_elements()->get<FiniteElement_name_mi_tag>().end()) {
312
313 for (auto peit = in_tets.pair_begin(); peit != in_tets.pair_end(); ++peit) {
314
315 auto lo =
316 prb_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
318 peit->first, (*fe_miit)->getFEUId()));
319 auto hi =
320 prb_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
322 peit->second, (*fe_miit)->getFEUId()));
323 numered_fes->insert(lo, hi);
324
326 std::cerr << "numered elements:" << std::endl;
327 for (; lo != hi; ++lo)
329 std::cerr << **lo << endl;
330 }
332 std::cerr << std::endl;
333
334 if (auto fe_ptr = data_ptr->feMethodPtr.lock()) {
335
337 "Number elements %d to evaluate at proc %d",
338 numered_fes->size(), m_field.get_comm_rank());
340
341 if (!cache_ptr) {
342 cache_ptr = boost::make_shared<CacheTuple>();
343 CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
344
345 MOFEM_LOG(
"FieldEvaluatorSync", Sev::noisy)
346 << "If you call that function many times in the loop consider to "
347 "set "
348 "cache_ptr outside of the loop. Otherwise code can be slow.";
349 }
350
351 CHKERR m_field.loop_finite_elements(prb_ptr, finite_element, *fe_ptr,
352 lower_rank, upper_rank, numered_fes,
353 bh, cache_ptr, verb);
354
355 } else
357 "Pointer to element does not exists");
358 }
359
361}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
@ MOFEM_DATA_INCONSISTENCY
FTensor::Index< 'n', SPACE_DIM > n
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< Ent_mi_tag >, 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.
UId getLocalUniqueIdCalculate() const
Generate UId for finite element entity.