MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
mesh_extras.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 "mesh_extras.hpp"
13 
14 using namespace std;
15 
16 namespace mfem
17 {
18 
19 namespace common
20 {
21 
22 ElementMeshStream::ElementMeshStream(Element::Type e)
23 {
24  *this << "MFEM mesh v1.0" << endl;
25  switch (e)
26  {
27  case Element::SEGMENT:
28  *this << "dimension" << endl << 1 << endl
29  << "elements" << endl << 1 << endl
30  << "1 1 0 1" << endl
31  << "boundary" << endl << 2 << endl
32  << "1 0 0" << endl
33  << "1 0 1" << endl
34  << "vertices" << endl
35  << 2 << endl
36  << 1 << endl
37  << 0 << endl
38  << 1 << endl;
39  break;
40  case Element::TRIANGLE:
41  *this << "dimension" << endl << 2 << endl
42  << "elements" << endl << 1 << endl
43  << "1 2 0 1 2" << endl
44  << "boundary" << endl << 3 << endl
45  << "1 1 0 1" << endl
46  << "1 1 1 2" << endl
47  << "1 1 2 0" << endl
48  << "vertices" << endl
49  << "3" << endl
50  << "2" << endl
51  << "0 0" << endl
52  << "1 0" << endl
53  << "0 1" << endl;
54  break;
55  case Element::QUADRILATERAL:
56  *this << "dimension" << endl << 2 << endl
57  << "elements" << endl << 1 << endl
58  << "1 3 0 1 2 3" << endl
59  << "boundary" << endl << 4 << endl
60  << "1 1 0 1" << endl
61  << "1 1 1 2" << endl
62  << "1 1 2 3" << endl
63  << "1 1 3 0" << endl
64  << "vertices" << endl
65  << "4" << endl
66  << "2" << endl
67  << "0 0" << endl
68  << "1 0" << endl
69  << "1 1" << endl
70  << "0 1" << endl;
71  break;
72  case Element::TETRAHEDRON:
73  *this << "dimension" << endl << 3 << endl
74  << "elements" << endl << 1 << endl
75  << "1 4 0 1 2 3" << endl
76  << "boundary" << endl << 4 << endl
77  << "1 2 0 2 1" << endl
78  << "1 2 1 2 3" << endl
79  << "1 2 2 0 3" << endl
80  << "1 2 0 1 3" << endl
81  << "vertices" << endl
82  << "4" << endl
83  << "3" << endl
84  << "0 0 0" << endl
85  << "1 0 0" << endl
86  << "0 1 0" << endl
87  << "0 0 1" << endl;
88  break;
89  case Element::HEXAHEDRON:
90  *this << "dimension" << endl << 3 << endl
91  << "elements" << endl << 1 << endl
92  << "1 5 0 1 2 3 4 5 6 7" << endl
93  << "boundary" << endl << 6 << endl
94  << "1 3 0 3 2 1" << endl
95  << "1 3 4 5 6 7" << endl
96  << "1 3 0 1 5 4" << endl
97  << "1 3 1 2 6 5" << endl
98  << "1 3 2 3 7 6" << endl
99  << "1 3 3 0 4 7" << endl
100  << "vertices" << endl
101  << "8" << endl
102  << "3" << endl
103  << "0 0 0" << endl
104  << "1 0 0" << endl
105  << "1 1 0" << endl
106  << "0 1 0" << endl
107  << "0 0 1" << endl
108  << "1 0 1" << endl
109  << "1 1 1" << endl
110  << "0 1 1" << endl;
111  break;
112  case Element::WEDGE:
113  *this << "dimension" << endl << 3 << endl
114  << "elements" << endl << 1 << endl
115  << "1 6 0 1 2 3 4 5" << endl
116  << "boundary" << endl << 5 << endl
117  << "1 2 2 1 0" << endl
118  << "1 2 3 4 5" << endl
119  << "1 3 0 1 4 3" << endl
120  << "1 3 1 2 5 4" << endl
121  << "1 3 2 0 3 5" << endl
122  << "vertices" << endl
123  << "6" << endl
124  << "3" << endl
125  << "0 0 0" << endl
126  << "1 0 0" << endl
127  << "0 1 0" << endl
128  << "0 0 1" << endl
129  << "1 0 1" << endl
130  << "0 1 1" << endl;
131  break;
132  default:
133  mfem_error("Invalid element type!");
134  break;
135  }
136 
137 }
138 
139 void
140 MergeMeshNodes(Mesh * mesh, int logging)
141 {
142  int dim = mesh->Dimension();
143  int sdim = mesh->SpaceDimension();
144 
145  double h_min, h_max, k_min, k_max;
146  mesh->GetCharacteristics(h_min, h_max, k_min, k_max);
147 
148  // Set tolerance for merging vertices
149  double tol = 1.0e-8 * h_min;
150 
151  if ( logging > 0 )
152  cout << "Euler Number of Initial Mesh: "
153  << ((dim==3)?mesh->EulerNumber() :
154  ((dim==2)?mesh->EulerNumber2D() :
155  mesh->GetNV() - mesh->GetNE())) << endl;
156 
157  vector<int> v2v(mesh->GetNV());
158 
159  Vector vd(sdim);
160 
161  for (int i = 0; i < mesh->GetNV(); i++)
162  {
163  Vector vi(mesh->GetVertex(i), sdim);
164 
165  v2v[i] = -1;
166 
167  for (int j = 0; j < i; j++)
168  {
169  Vector vj(mesh->GetVertex(j), sdim);
170  add(vi, -1.0, vj, vd);
171 
172  if ( vd.Norml2() < tol )
173  {
174  v2v[i] = j;
175  break;
176  }
177  }
178  if ( v2v[i] < 0 ) { v2v[i] = i; }
179  }
180 
181  // renumber elements
182  for (int i = 0; i < mesh->GetNE(); i++)
183  {
184  Element *el = mesh->GetElement(i);
185  int *v = el->GetVertices();
186  int nv = el->GetNVertices();
187  for (int j = 0; j < nv; j++)
188  {
189  v[j] = v2v[v[j]];
190  }
191  }
192  // renumber boundary elements
193  for (int i = 0; i < mesh->GetNBE(); i++)
194  {
195  Element *el = mesh->GetBdrElement(i);
196  int *v = el->GetVertices();
197  int nv = el->GetNVertices();
198  for (int j = 0; j < nv; j++)
199  {
200  v[j] = v2v[v[j]];
201  }
202  }
203 
204  mesh->RemoveUnusedVertices();
205 
206  if ( logging > 0 )
207  {
208  cout << "Euler Number of Final Mesh: "
209  << ((dim==3) ? mesh->EulerNumber() :
210  ((dim==2) ? mesh->EulerNumber2D() :
211  mesh->GetNV() - mesh->GetNE()))
212  << endl;
213  }
214 }
215 
216 void AttrToMarker(int max_attr, const Array<int> &attrs, Array<int> &marker)
217 {
218  MFEM_ASSERT(attrs.Max() <= max_attr, "Invalid attribute number present.");
219 
220  marker.SetSize(max_attr);
221  if (attrs.Size() == 1 && attrs[0] == -1)
222  {
223  marker = 1;
224  }
225  else
226  {
227  marker = 0;
228  for (int j=0; j<attrs.Size(); j++)
229  {
230  int attr = attrs[j];
231  MFEM_VERIFY(attr > 0, "Attribute number less than one!");
232  marker[attr-1] = 1;
233  }
234  }
235 }
236 
237 void KershawTransformation::Eval(Vector &V, ElementTransformation &T,
238  const IntegrationPoint &ip)
239 {
240  V = 0.0;
241  Vector pos(dim);
242  T.Transform(ip, pos);
243  double x = pos(0), y = pos(1), z = dim == 3 ? pos(2) : 0;
244  double X, Y, Z;
245 
246  X = x;
247 
248  int layer = x*6.0;
249  double lambda = (x-layer/6.0)*6;
250 
251  // The x-range is split in 6 layers going from left-to-left, left-to-right,
252  // right-to-left (2 layers), left-to-right and right-to-right yz-faces.
253  switch (layer)
254  {
255  case 0:
256  Y = left(epsy, y);
257  Z = left(epsz, z);
258  break;
259  case 1:
260  case 4:
261  Y = step(left(epsy, y), right(epsy, y), lambda);
262  Z = step(left(epsz, z), right(epsz, z), lambda);
263  break;
264  case 2:
265  Y = step(right(epsy, y), left(epsy, y), lambda/2);
266  Z = step(right(epsz, z), left(epsz, z), lambda/2);
267  break;
268  case 3:
269  Y = step(right(epsy, y), left(epsy, y), (1+lambda)/2);
270  Z = step(right(epsz, z), left(epsz, z), (1+lambda)/2);
271  break;
272  default:
273  Y = right(epsy, y);
274  Z = right(epsz, z);
275  break;
276  }
277  V.SetSize(dim);
278  V(0) = X;
279  V(1) = Y;
280  if (dim == 3) { V(2) = Z; }
281 }
282 
283 } // namespace common
284 
285 } // namespace mfem
mfem::Array::Size
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
mfem::Mesh::GetVertex
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1012
mfem::Element::GetVertices
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
mfem::common::AttrToMarker
void AttrToMarker(int max_attr, const Array< int > &attrs, Array< int > &marker)
Convert a set of attribute numbers to a marker array.
Definition: mesh_extras.cpp:216
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
mfem::Mesh::GetNBE
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:926
mfem::Mesh::GetCharacteristics
void GetCharacteristics(double &h_min, double &h_max, double &kappa_min, double &kappa_max, Vector *Vh=NULL, Vector *Vk=NULL)
Definition: mesh.cpp:193
mfem::Mesh::EulerNumber2D
int EulerNumber2D() const
Equals 1 - num_holes.
Definition: mesh.hpp:1003
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:923
mfem::common::MergeMeshNodes
void MergeMeshNodes(Mesh *mesh, int logging)
Merges vertices which lie at the same location.
Definition: mesh_extras.cpp:140
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:313
mesh_extras.hpp
mfem::mfem_error
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
mfem::Element::Type
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
mfem::Array::Max
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
mfem::Mesh::GetElement
const Element * GetElement(int i) const
Definition: mesh.hpp:1040
mfem::Mesh::Dimension
int Dimension() const
Definition: mesh.hpp:1006
mfem::Mesh::SpaceDimension
int SpaceDimension() const
Definition: mesh.hpp:1007
mfem::Mesh::RemoveUnusedVertices
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:11678
mfem::Array::SetSize
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
mfem::Mesh::EulerNumber
int EulerNumber() const
Equals 1 + num_holes - num_loops.
Definition: mesh.hpp:1000
mfem::Mesh::GetNV
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:920
mfem::IntegrationPoint
Class for integration point with weight.
Definition: intrules.hpp:25
dim
int dim
Definition: ex24.cpp:53
mfem::Element::GetNVertices
virtual int GetNVertices() const =0
mfem::Vector
Vector data type.
Definition: vector.hpp:60
mfem::ElementTransformation::Transform
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
mfem::Element
Abstract data type element.
Definition: element.hpp:28
mfem::Mesh::GetBdrElement
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1044
Generated on Wed Dec 7 2022 10:41:49 for MFEM by   doxygen 1.8.5