| 336 | |
| 337 | template<int DIM> |
| 338 | void T_GapFunction<DIM> :: Evaluate(const BaseMappedIntegrationPoint & ip, |
| 339 | FlatVector<> result) const |
| 340 | { |
| 341 | LocalHeapMem<100000> lh("gapfunction"); |
| 342 | auto & trafo1 = ip.GetTransformation(); |
| 343 | const auto & el1 = ma->GetElement(trafo1.GetElementId()); |
| 344 | result = 0; |
| 345 | if (!master.Mask().Test(el1.GetIndex())) return; |
| 346 | |
| 347 | // int intorder2 = 2*displacement->GetFESpace()->GetOrder(); |
| 348 | |
| 349 | auto & trafo1_def = trafo1.AddDeformation(displacement.get(), lh); |
| 350 | |
| 351 | double inv_fac = GetDomIn(*ma, el1) == 0 ? -1. : 1.; |
| 352 | |
| 353 | auto & ip1 = ip.IP(); |
| 354 | Vec<DIM> p1; |
| 355 | trafo1_def.CalcPoint(ip1, p1); |
| 356 | |
| 357 | double mindist = 1e99; |
| 358 | // result = std::numeric_limits<double>::infinity(); |
| 359 | result = std::numeric_limits<double>::max(); |
| 360 | |
| 361 | // find all bound-2 elements closer to p1 than h |
| 362 | netgen::Point<DIM> ngp1; |
| 363 | for (int j = 0; j < DIM; j++) |
| 364 | ngp1(j) = p1(j); |
| 365 | |
| 366 | auto & mip = static_cast<const DimMappedIntegrationPoint<DIM>&>(ip); |
| 367 | |
| 368 | auto hcurrent = h/(1024.*1024.); |
| 369 | int found = 2; |
| 370 | while(found>0 && hcurrent<=h) |
| 371 | { |
| 372 | netgen::Box<DIM> box(ngp1, ngp1); |
| 373 | box.Increase(hcurrent); |
| 374 | |
| 375 | searchtree->GetFirstIntersecting |
| 376 | (box.PMin(), box.PMax(), |
| 377 | [&] (int elnr2) |
| 378 | { |
| 379 | auto el2 = ma->GetElement( ElementId (BND, elnr2) ); |
| 380 | double inv_fac2 = GetDomIn(*ma, el2) == 0 ? -1. : 1.; |
| 381 | HeapReset hr(lh); |
| 382 | |
| 383 | bool common_vertex = false; |
| 384 | for (auto s_v : el1.Vertices() ) |
| 385 | for (auto v : el2.Vertices() ) |
| 386 | if(s_v==v) |
| 387 | common_vertex = true; |
| 388 | if (common_vertex) return false; |
| 389 | auto & trafo2 = ma->GetTrafo (el2, lh); |
| 390 | auto & trafo2_def = trafo2.AddDeformation(displacement.get(), lh); |
| 391 | |
| 392 | Vec<DIM> p2; |
| 393 | IntegrationPoint ip2; |
| 394 | double dist = FindClosestPoint<DIM-1,DIM>(p1, inv_fac * inv_fac2 * mip.GetNV(), mindist, trafo2_def, ip2, p2, both_sides); |
| 395 | if(dist<mindist && dist < h) |
nothing calls this directly
no test coverage detected