MFEM  v4.6.0
Finite element discretization library
kdtree.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 
14 namespace mfem
15 {
16 
17 template<>
18 void KDTreeNodalProjection<2>::Project(const Vector& coords,const Vector& src,
19  int ordering, double 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  double 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 
84 template<>
85 void KDTreeNodalProjection<3>::Project(const Vector& coords,const Vector& src,
86  int ordering, double 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  double 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 
152 template<>
153 void KDTreeNodalProjection<2>::Project(const GridFunction& gf, double 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  double 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 
277 template<>
278 void KDTreeNodalProjection<3>::Project(const GridFunction& gf, double 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  double 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::FiniteElement
Abstract class for all finite elements.
Definition: fe_base.hpp:233
mfem::KDTree::PointND::xx
Tfloat xx[ndim]
Coordinates of the point.
Definition: kdtree.hpp:112
trans
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
mfem::GridFunction
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
mfem::FiniteElementSpace::GetFE
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
mfem::GridFunction::FESpace
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
mfem::FiniteElementSpace::GetNE
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
mfem::FiniteElementSpace::GetVDim
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
mfem::FiniteElementSpace::GetElementVDofs
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
p
double p(const Vector &x, double t)
Definition: navier_mms.cpp:53
mfem::FiniteElementSpace::GetOrdering
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:721
dim
int dim
Definition: ex24.cpp:53
mfem::FiniteElementSpace::GetVSize
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
mfem::Vector
Vector data type.
Definition: vector.hpp:58
mfem::DenseMatrix::SetSize
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
mfem::FiniteElement::GetNodes
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:389
mfem::KDTreeNodalProjection::Project
virtual void Project(const Vector &coords, const Vector &src, int ordering=Ordering::byNODES, double lerr=1e-8)
mfem::FiniteElementSpace::GetElementTransformation
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:769
Generated on Wed Sep 27 2023 14:52:30 for MFEM by   doxygen 1.8.14