| 82 | |
| 83 | template<int DIMS, int DIMR> |
| 84 | double FindClosestPoint( Vec<DIMR> pmaster, Vec<DIMR> n, double h, const ElementTransformation & trafo, IntegrationPoint & ip, Vec<DIMR> & p, bool both_sides) |
| 85 | // input arguments: pmaster, n, h (maximum distance), trafo |
| 86 | // output arguments: ip, p |
| 87 | { |
| 88 | Vec<DIMS> min_lam = 0; |
| 89 | double min_dist = 1e99; |
| 90 | |
| 91 | min_lam = 1./(DIMS+1); |
| 92 | |
| 93 | // Todo: line search, stop criterion |
| 94 | for([[maybe_unused]] auto i : Range(4) ) |
| 95 | { |
| 96 | ip = min_lam; |
| 97 | MappedIntegrationPoint<DIMS, DIMR> mip{ip, trafo}; |
| 98 | T2<DIMS> t2{mip, pmaster}; |
| 99 | bool is_front = InnerProduct(n, mip.GetNV()) < 0; |
| 100 | if(both_sides) |
| 101 | is_front = true; |
| 102 | |
| 103 | if constexpr (DIMS==1) |
| 104 | { |
| 105 | // check end points |
| 106 | for(double lam : {0.,1.}) |
| 107 | { |
| 108 | auto dist = t2(lam); |
| 109 | if(is_front && dist<min_dist) |
| 110 | { |
| 111 | min_dist = dist; |
| 112 | min_lam = lam; |
| 113 | } |
| 114 | } |
| 115 | |
| 116 | auto lam = t2.CalcMinimum(); |
| 117 | auto dist = t2(lam); |
| 118 | if(is_front && lam[0]>0 && lam[0] < 1) |
| 119 | { |
| 120 | min_dist = dist; |
| 121 | min_lam = lam; |
| 122 | } |
| 123 | } |
| 124 | if constexpr (DIMS==2) |
| 125 | { |
| 126 | auto getDist = [&] ( auto l ) |
| 127 | { |
| 128 | Vec<DIMR> p; |
| 129 | trafo.CalcPoint( {l}, p ); |
| 130 | return L2Norm2( p-pmaster ); |
| 131 | }; |
| 132 | |
| 133 | // check corner points and edges separately |
| 134 | ArrayMem<Vec<DIMS>, 4> points; |
| 135 | auto eltype = trafo.GetElementType(); |
| 136 | if(eltype==ET_TRIG) |
| 137 | points = { {0,0}, {0,1}, {1,0} }; |
| 138 | else if(eltype==ET_QUAD) |
| 139 | points = { {0,0}, {0,1}, {1,1}, {1,0} }; |
| 140 | |
| 141 | for (Vec<DIMS> lam : points) |
nothing calls this directly
no test coverage detected