MFEM  v4.4.0
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 } // namespace common
238 
239 } // namespace mfem
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1005
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void AttrToMarker(int max_attr, const Array< int > &attrs, Array< int > &marker)
Convert a set of attribute numbers to a marker array.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:931
void GetCharacteristics(double &h_min, double &h_max, double &kappa_min, double &kappa_max, Vector *Vh=NULL, Vector *Vk=NULL)
Definition: mesh.cpp:193
int EulerNumber2D() const
Equals 1 - num_holes.
Definition: mesh.hpp:996
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
void MergeMeshNodes(Mesh *mesh, int logging)
Merges vertices which lie at the same location.
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:300
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
Type
Constants for the classes derived from Element.
Definition: element.hpp:41
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
const Element * GetElement(int i) const
Definition: mesh.hpp:1033
int Dimension() const
Definition: mesh.hpp:999
int SpaceDimension() const
Definition: mesh.hpp:1000
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:11600
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
int EulerNumber() const
Equals 1 + num_holes - num_loops.
Definition: mesh.hpp:993
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:925
int dim
Definition: ex24.cpp:53
virtual int GetNVertices() const =0
Vector data type.
Definition: vector.hpp:60
Abstract data type element.
Definition: element.hpp:28
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1037