MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
kdtree.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12#include "kdtree.hpp"
13
14namespace mfem
15{
16
17template<>
18void KDTreeNodalProjection<2>::Project(const Vector& coords,const Vector& src,
19 int ordering, real_t lerr)
20{
21 const int dim=dest->FESpace()->GetMesh()->SpaceDimension();
22 const int vd=dest->VectorDim(); // dimension of the vector field
23 const int np=src.Size()/vd; // number of points
24 int ind;
25 real_t dist;
26 bool pt_inside_bbox;
27 KDTree2D::PointND pnd;
28 for (int i=0; i<np; i++)
29 {
30 pnd.xx[0]=coords(i*dim+0);
31 pnd.xx[1]=coords(i*dim+1);
32
33 pt_inside_bbox=true;
34 for (int j=0; j<dim; j++)
35 {
36 if (pnd.xx[j]>(maxbb[j]+lerr)) {pt_inside_bbox=false; break;}
37 if (pnd.xx[j]<(minbb[j]-lerr)) {pt_inside_bbox=false; break;}
38 }
39
40 if (pt_inside_bbox)
41 {
42 kdt->FindClosestPoint(pnd,ind,dist);
43 if (dist<lerr)
44 {
45 if (dest->FESpace()->GetOrdering()==Ordering::byNODES)
46 {
47 if (ordering==Ordering::byNODES)
48 {
49 for (int di=0; di<vd; di++)
50 {
51 (*dest)[di*np+ind]=src[di*np+i];
52 }
53 }
54 else
55 {
56 for (int di=0; di<vd; di++)
57 {
58 (*dest)[di*np+ind]=src[di+i*vd];
59 }
60 }
61 }
62 else
63 {
64 if (ordering==Ordering::byNODES)
65 {
66 for (int di=0; di<vd; di++)
67 {
68 (*dest)[di+ind*vd]=src[di*np+i];
69 }
70 }
71 else
72 {
73 for (int di=0; di<vd; di++)
74 {
75 (*dest)[di+ind*vd]=src[di+i*vd];
76 }
77 }
78 }
79 }
80 }
81 }
82}
83
84template<>
85void KDTreeNodalProjection<3>::Project(const Vector& coords,const Vector& src,
86 int ordering, real_t lerr)
87{
88 const int dim=dest->FESpace()->GetMesh()->SpaceDimension();
89 const int vd=dest->VectorDim(); // dimension of the vector field
90 const int np=src.Size()/vd; // number of points
91 int ind;
92 real_t dist;
93 bool pt_inside_bbox;
94 KDTree3D::PointND pnd;
95 for (int i=0; i<np; i++)
96 {
97 pnd.xx[0]=coords(i*dim+0);
98 pnd.xx[1]=coords(i*dim+1);
99 pnd.xx[2]=coords(i*dim+2);
100
101 pt_inside_bbox=true;
102 for (int j=0; j<dim; j++)
103 {
104 if (pnd.xx[j]>(maxbb[j]+lerr)) {pt_inside_bbox=false; break;}
105 if (pnd.xx[j]<(minbb[j]-lerr)) {pt_inside_bbox=false; break;}
106 }
107
108 if (pt_inside_bbox)
109 {
110 kdt->FindClosestPoint(pnd,ind,dist);
111 if (dist<lerr)
112 {
113 if (dest->FESpace()->GetOrdering()==Ordering::byNODES)
114 {
115 if (ordering==Ordering::byNODES)
116 {
117 for (int di=0; di<vd; di++)
118 {
119 (*dest)[di*np+ind]=src[di*np+i];
120 }
121 }
122 else
123 {
124 for (int di=0; di<vd; di++)
125 {
126 (*dest)[di*np+ind]=src[di+i*vd];
127 }
128 }
129 }
130 else
131 {
132 if (ordering==Ordering::byNODES)
133 {
134 for (int di=0; di<vd; di++)
135 {
136 (*dest)[di+ind*vd]=src[di*np+i];
137 }
138 }
139 else
140 {
141 for (int di=0; di<vd; di++)
142 {
143 (*dest)[di+ind*vd]=src[di+i*vd];
144 }
145 }
146 }
147 }
148 }
149 }
150}
151
152template<>
153void KDTreeNodalProjection<2>::Project(const GridFunction& gf, real_t lerr)
154{
155 int ordering = gf.FESpace()->GetOrdering();
156 Vector coo;
157 int np=gf.FESpace()->GetVSize()/gf.FESpace()->GetVDim();
158 coo.SetSize(np*2);
159 int vd=dest->VectorDim();
160 int ind;
161 real_t dist;
162
163 Vector maxbb_src(2);
164 Vector minbb_src(2);
165
166 // extract the nodal coordinates from gf
167 {
168 ElementTransformation *trans;
169 const IntegrationRule* ir=nullptr;
170 Array<int> vdofs;
171 DenseMatrix elco;
172 int isca=1;
173 if (gf.FESpace()->GetOrdering()==Ordering::byVDIM)
174 {
175 isca=gf.FESpace()->GetVDim();
176 }
177
178 // initialize bbmax and bbmin
179 const FiniteElement* el=gf.FESpace()->GetFE(0);
180 trans = gf.FESpace()->GetElementTransformation(0);
181 ir=&(el->GetNodes());
182 gf.FESpace()->GetElementVDofs(0,vdofs);
183 elco.SetSize(2,ir->GetNPoints());
184 trans->Transform(*ir,elco);
185 for (int d=0; d<2; d++)
186 {
187 maxbb_src(d)=elco(d,0);
188 minbb_src(d)=elco(d,0);
189 }
190
191 for (int i=0; i<gf.FESpace()->GetNE(); i++)
192 {
193 el=gf.FESpace()->GetFE(i);
194 //get the element transformation
195 trans = gf.FESpace()->GetElementTransformation(i);
196 ir=&(el->GetNodes());
197 gf.FESpace()->GetElementVDofs(i,vdofs);
198 elco.SetSize(2,ir->GetNPoints());
199 trans->Transform(*ir,elco);
200 for (int p=0; p<ir->GetNPoints(); p++)
201 {
202 for (int d=0; d<2; d++)
203 {
204 coo[vdofs[p]*2/isca+d]=elco(d,p);
205
206 if (maxbb_src(d)<elco(d,p)) {maxbb_src(d)=elco(d,p);}
207 if (minbb_src(d)>elco(d,p)) {minbb_src(d)=elco(d,p);}
208 }
209 }
210 }
211 }
212
213 maxbb_src+=lerr;
214 minbb_src-=lerr;
215
216 // check for intersection
217 bool flag;
218 {
219 flag=true;
220 for (int i=0; i<2; i++)
221 {
222 if (minbb_src(i)>maxbb(i)) {flag=false;}
223 if (maxbb_src(i)<minbb(i)) {flag=false;}
224 }
225 if (flag==false) {return;}
226 }
227
228 {
229 KDTree2D::PointND pnd;
230 for (int i=0; i<np; i++)
231 {
232 pnd.xx[0]=coo(i*2+0);
233 pnd.xx[1]=coo(i*2+1);
234
235 kdt->FindClosestPoint(pnd,ind,dist);
236 if (dist<lerr)
237 {
238 if (dest->FESpace()->GetOrdering()==Ordering::byNODES)
239 {
240 if (ordering==Ordering::byNODES)
241 {
242 for (int di=0; di<vd; di++)
243 {
244 (*dest)[di*np+ind]=gf[di*np+i];
245 }
246 }
247 else
248 {
249 for (int di=0; di<vd; di++)
250 {
251 (*dest)[di*np+ind]=gf[di+i*vd];
252 }
253 }
254 }
255 else
256 {
257 if (ordering==Ordering::byNODES)
258 {
259 for (int di=0; di<vd; di++)
260 {
261 (*dest)[di+ind*vd]=gf[di*np+i];
262 }
263 }
264 else
265 {
266 for (int di=0; di<vd; di++)
267 {
268 (*dest)[di+ind*vd]=gf[di+i*vd];
269 }
270 }
271 }
272 }
273 }
274 }
275}
276
277template<>
278void KDTreeNodalProjection<3>::Project(const GridFunction& gf, real_t lerr)
279{
280 int ordering = gf.FESpace()->GetOrdering();
281 int dim=dest->FESpace()->GetMesh()->SpaceDimension();
282 Vector coo;
283 int np=gf.FESpace()->GetVSize()/gf.FESpace()->GetVDim();
284 coo.SetSize(np*dim);
285 int vd=dest->VectorDim();
286 int ind;
287 real_t dist;
288
289 Vector maxbb_src(dim);
290 Vector minbb_src(dim);
291
292 // extract the nodal coordinates from gf
293 {
294 ElementTransformation *trans;
295 const IntegrationRule* ir=nullptr;
296 Array<int> vdofs;
297 DenseMatrix elco;
298 int isca=1;
299 if (gf.FESpace()->GetOrdering()==Ordering::byVDIM)
300 {
301 isca=gf.FESpace()->GetVDim();
302 }
303
304 // initialize bbmax and bbmin
305 const FiniteElement* el=gf.FESpace()->GetFE(0);
306 trans = gf.FESpace()->GetElementTransformation(0);
307 ir=&(el->GetNodes());
308 gf.FESpace()->GetElementVDofs(0,vdofs);
309 elco.SetSize(dim,ir->GetNPoints());
310 trans->Transform(*ir,elco);
311 for (int d=0; d<dim; d++)
312 {
313 maxbb_src(d)=elco(d,0);
314 minbb_src(d)=elco(d,0);
315 }
316
317 for (int i=0; i<gf.FESpace()->GetNE(); i++)
318 {
319 el=gf.FESpace()->GetFE(i);
320 // get the element transformation
321 trans = gf.FESpace()->GetElementTransformation(i);
322 ir=&(el->GetNodes());
323 gf.FESpace()->GetElementVDofs(i,vdofs);
324 elco.SetSize(dim,ir->GetNPoints());
325 trans->Transform(*ir,elco);
326 for (int p=0; p<ir->GetNPoints(); p++)
327 {
328 for (int d=0; d<dim; d++)
329 {
330 coo[vdofs[p]*dim/isca+d]=elco(d,p);
331
332 if (maxbb_src(d)<elco(d,p)) {maxbb_src(d)=elco(d,p);}
333 if (minbb_src(d)>elco(d,p)) {minbb_src(d)=elco(d,p);}
334 }
335 }
336 }
337 }
338
339 maxbb_src+=lerr;
340 minbb_src-=lerr;
341
342 // check for intersection
343 bool flag;
344 {
345 flag=true;
346 for (int i=0; i<dim; i++)
347 {
348 if (minbb_src(i)>maxbb(i)) {flag=false;}
349 if (maxbb_src(i)<minbb(i)) {flag=false;}
350 }
351 if (flag==false) {return;}
352 }
353
354 {
355 KDTree3D::PointND pnd;
356 for (int i=0; i<np; i++)
357 {
358 pnd.xx[0]=coo(i*dim+0);
359 pnd.xx[1]=coo(i*dim+1);
360 pnd.xx[2]=coo(i*dim+2);
361
362 kdt->FindClosestPoint(pnd,ind,dist);
363 if (dist<lerr)
364 {
365 if (dest->FESpace()->GetOrdering()==Ordering::byNODES)
366 {
367 if (ordering==Ordering::byNODES)
368 {
369 for (int di=0; di<vd; di++)
370 {
371 (*dest)[di*np+ind]=gf[di*np+i];
372 }
373 }
374 else
375 {
376 for (int di=0; di<vd; di++)
377 {
378 (*dest)[di*np+ind]=gf[di+i*vd];
379 }
380 }
381 }
382 else
383 {
384 if (ordering==Ordering::byNODES)
385 {
386 for (int di=0; di<vd; di++)
387 {
388 (*dest)[di+ind*vd]=gf[di*np+i];
389 }
390 }
391 else
392 {
393 for (int di=0; di<vd; di++)
394 {
395 (*dest)[di+ind*vd]=gf[di+i*vd];
396 }
397 }
398 }
399 }
400 }
401 }
402}
403
404} // namespace mfem
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
mfem::FiniteElementSpace::GetElementTransformation
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition fespace.hpp:773
mfem::FiniteElementSpace::GetElementVDofs
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition fespace.cpp:287
mfem::FiniteElementSpace::GetFE
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3168
mfem::FiniteElementSpace::GetOrdering
Ordering::Type GetOrdering() const
Return the ordering method.
Definition fespace.hpp:725
mfem::FiniteElementSpace::GetVSize
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
mfem::FiniteElement
Abstract class for all finite elements.
Definition fe_base.hpp:239
mfem::FiniteElement::GetNodes
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:395
mfem::GridFunction
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
mfem::GridFunction::FESpace
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
mfem::KDTreeNodalProjection::Project
virtual void Project(const Vector &coords, const Vector &src, int ordering=Ordering::byNODES, real_t lerr=1e-8)
mfem::Vector
Vector data type.
Definition vector.hpp:80
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
dim
int dim
Definition ex24.cpp:53
trans
void trans(const Vector &u, Vector &x)
Definition ex27.cpp:412
mfem::real_t
float real_t
Definition config.hpp:43
p
real_t p(const Vector &x, real_t t)
Definition navier_mms.cpp:53
mfem::KDTree::PointND::xx
Tfloat xx[ndim]
Coordinates of the point.
Definition kdtree.hpp:133
doxygen 1.11.0