179 {
180 CoreInterface &m_field =
cOre;
182
183 std::vector<EntityHandle> leafs_out;
184 CHKERR data_ptr->treePtr->distance_search(point, distance, leafs_out);
186 for (auto lit : leafs_out)
187 CHKERR m_field.get_moab().get_entities_by_dimension(lit,
D, tree_ents,
188 true);
189
191 MOFEM_LOG(
"SELF", Sev::noisy) <<
"tree entities: " << tree_ents;
192
193 data_ptr->evalPointEntityHandle.resize(data_ptr->nbEvalPoints);
194 std::fill(data_ptr->evalPointEntityHandle.begin(),
195 data_ptr->evalPointEntityHandle.end(), 0);
196
197 std::vector<double> local_coords;
198 std::vector<double> shape;
199
200 auto nb_eval_points = data_ptr->nbEvalPoints;
201
202 for (auto tet : tree_ents) {
206 "Wrong element type, Field Evaluator not implemented for Quads and Hexes");
207 }
209 int num_nodes;
210 CHKERR m_field.get_moab().get_connectivity(tet, conn, num_nodes,
true);
211
212 if constexpr (
D == 3) {
213
214 local_coords.resize(3 * nb_eval_points);
215 shape.resize(4 * nb_eval_points);
216
217 std::array<double, 12> coords;
218 CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
219
221 coords.data(), &data_ptr->evalPoints[0], nb_eval_points,
222 &local_coords[0]);
224 &local_coords[1], &local_coords[2],
225 nb_eval_points);
226
229 &shape[0], &shape[1], &shape[2], &shape[3]};
231 &data_ptr->shapeFunctions(0, 0), &data_ptr->shapeFunctions(0, 1),
232 &data_ptr->shapeFunctions(0, 2), &data_ptr->shapeFunctions(0, 3)};
233
236 &local_coords[0], &local_coords[1], &local_coords[2]};
238 &(data_ptr->localCoords(0, 0)), &(data_ptr->localCoords(0, 1)),
239 &(data_ptr->localCoords(0, 2))};
240
241 for (
int n = 0;
n != nb_eval_points; ++
n) {
242
243 const double eps = data_ptr->eps;
244 if (t_shape(0) >= 0 -
eps && t_shape(0) <= 1 +
eps &&
245
246 t_shape(1) >= 0 -
eps && t_shape(1) <= 1 +
eps &&
247
248 t_shape(2) >= 0 -
eps && t_shape(2) <= 1 +
eps &&
249
250 t_shape(3) >= 0 -
eps && t_shape(3) <= 1 +
eps) {
251
252 data_ptr->evalPointEntityHandle[
n] = tet;
253 t_shape_data(i4) = t_shape(i4);
254 t_local_data(j3) = t_local(j3);
255 }
256
257 ++t_shape;
258 ++t_shape_data;
259 ++t_local;
260 ++t_local_data;
261 }
262 }
263
264 if constexpr (
D == 2) {
265
266 local_coords.resize(2 * nb_eval_points);
267 shape.resize(3 * nb_eval_points);
268
269 std::array<double, 9> coords;
270 CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords.data());
271
273 coords.data(), &data_ptr->evalPoints[0], nb_eval_points,
274 &local_coords[0]);
276 &local_coords[1], nb_eval_points);
277
280 &shape[0], &shape[1], &shape[2]};
282 &data_ptr->shapeFunctions(0, 0), &data_ptr->shapeFunctions(0, 1),
283 &data_ptr->shapeFunctions(0, 2)};
284
287 &local_coords[0], &local_coords[1]};
288 data_ptr->localCoords.resize(nb_eval_points, 2);
290 &(data_ptr->localCoords(0, 0)), &(data_ptr->localCoords(0, 1))};
291
292 for (
int n = 0;
n != nb_eval_points; ++
n) {
293
294 const double eps = data_ptr->eps;
295 if (t_shape(0) >= 0 -
eps && t_shape(0) <= 1 +
eps &&
296
297 t_shape(1) >= 0 -
eps && t_shape(1) <= 1 +
eps &&
298
299 t_shape(2) >= 0 -
eps && t_shape(2) <= 1 +
eps) {
300
301 data_ptr->evalPointEntityHandle[
n] = tet;
302 t_shape_data(i3) = t_shape(i3);
303 t_local_data(j2) = t_local(j2);
304 }
305
306 ++t_shape;
307 ++t_shape_data;
308 ++t_local;
309 ++t_local_data;
310 }
311 }
312 }
313
314 const Problem *prb_ptr;
315 CHKERR m_field.get_problem(problem, &prb_ptr);
316 boost::shared_ptr<NumeredEntFiniteElement_multiIndex> numered_fes(
318
320 in_tets.insert_list(data_ptr->evalPointEntityHandle.begin(),
321 data_ptr->evalPointEntityHandle.end());
322 in_tets = in_tets.subset_by_dimension(
D);
323
325 MOFEM_LOG(
"SELF", Sev::noisy) <<
"in tets: " << in_tets << endl;
326
327 auto fe_miit =
328 m_field.get_finite_elements()->get<FiniteElement_name_mi_tag>().find(
329 finite_element);
330 if (fe_miit !=
331 m_field.get_finite_elements()->get<FiniteElement_name_mi_tag>().end()) {
332
333 for (auto peit = in_tets.pair_begin(); peit != in_tets.pair_end(); ++peit) {
334
335 auto lo =
336 prb_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().lower_bound(
338 peit->first, (*fe_miit)->getFEUId()));
339 auto hi =
340 prb_ptr->numeredFiniteElementsPtr->get<Unique_mi_tag>().upper_bound(
342 peit->second, (*fe_miit)->getFEUId()));
343 numered_fes->insert(lo, hi);
344
346 std::cerr << "numered elements:" << std::endl;
347 for (; lo != hi; ++lo)
349 std::cerr << **lo << endl;
350 }
352 std::cerr << std::endl;
353
354 if (auto fe_ptr = data_ptr->feMethodPtr.lock()) {
355
357 "Number elements %d to evaluate at proc %d",
358 numered_fes->size(), m_field.get_comm_rank());
360
361 if (!cache_ptr) {
362 cache_ptr = boost::make_shared<CacheTuple>();
363 CHKERR m_field.cache_problem_entities(prb_ptr->getName(), cache_ptr);
364
365 MOFEM_LOG(
"FieldEvaluatorSync", Sev::noisy)
366 << "If you call that function many times in the loop consider to "
367 "set "
368 "cache_ptr outside of the loop. Otherwise code can be slow.";
369 }
370
371 CHKERR m_field.loop_finite_elements(prb_ptr, finite_element, *fe_ptr,
372 lower_rank, upper_rank, numered_fes,
373 bh, cache_ptr, verb);
374
375 } else
377 "Pointer to element does not exists");
378 }
379
381}
#define MOFEM_LOG_SYNCHRONISE(comm)
Synchronise "SYNC" channel.
#define MOFEM_LOG_C(channel, severity, format,...)
@ MOFEM_DATA_INCONSISTENCY
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.
const double n
refractive index of diffusive medium
auto type_from_handle(const EntityHandle h)
get type from entity handle
UId getLocalUniqueIdCalculate() const
Generate UId for finite element entity.
Public Member Functions inherited from MoFEM::UnknownInterface