MFEM  v4.3.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-2021, 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  default:
113  mfem_error("Invalid element type!");
114  break;
115  }
116 
117 }
118 
119 void
120 MergeMeshNodes(Mesh * mesh, int logging)
121 {
122  int dim = mesh->Dimension();
123  int sdim = mesh->SpaceDimension();
124 
125  double h_min, h_max, k_min, k_max;
126  mesh->GetCharacteristics(h_min, h_max, k_min, k_max);
127 
128  // Set tolerance for merging vertices
129  double tol = 1.0e-8 * h_min;
130 
131  if ( logging > 0 )
132  cout << "Euler Number of Initial Mesh: "
133  << ((dim==3)?mesh->EulerNumber() :
134  ((dim==2)?mesh->EulerNumber2D() :
135  mesh->GetNV() - mesh->GetNE())) << endl;
136 
137  vector<int> v2v(mesh->GetNV());
138 
139  Vector vd(sdim);
140 
141  for (int i = 0; i < mesh->GetNV(); i++)
142  {
143  Vector vi(mesh->GetVertex(i), sdim);
144 
145  v2v[i] = -1;
146 
147  for (int j = 0; j < i; j++)
148  {
149  Vector vj(mesh->GetVertex(j), sdim);
150  add(vi, -1.0, vj, vd);
151 
152  if ( vd.Norml2() < tol )
153  {
154  v2v[i] = j;
155  break;
156  }
157  }
158  if ( v2v[i] < 0 ) { v2v[i] = i; }
159  }
160 
161  // renumber elements
162  for (int i = 0; i < mesh->GetNE(); i++)
163  {
164  Element *el = mesh->GetElement(i);
165  int *v = el->GetVertices();
166  int nv = el->GetNVertices();
167  for (int j = 0; j < nv; j++)
168  {
169  v[j] = v2v[v[j]];
170  }
171  }
172  // renumber boundary elements
173  for (int i = 0; i < mesh->GetNBE(); i++)
174  {
175  Element *el = mesh->GetBdrElement(i);
176  int *v = el->GetVertices();
177  int nv = el->GetNVertices();
178  for (int j = 0; j < nv; j++)
179  {
180  v[j] = v2v[v[j]];
181  }
182  }
183 
184  mesh->RemoveUnusedVertices();
185 
186  if ( logging > 0 )
187  {
188  cout << "Euler Number of Final Mesh: "
189  << ((dim==3) ? mesh->EulerNumber() :
190  ((dim==2) ? mesh->EulerNumber2D() :
191  mesh->GetNV() - mesh->GetNE()))
192  << endl;
193  }
194 }
195 
196 void AttrToMarker(int max_attr, const Array<int> &attrs, Array<int> &marker)
197 {
198  MFEM_ASSERT(attrs.Max() <= max_attr, "Invalid attribute number present.");
199 
200  marker.SetSize(max_attr);
201  if (attrs.Size() == 1 && attrs[0] == -1)
202  {
203  marker = 1;
204  }
205  else
206  {
207  marker = 0;
208  for (int j=0; j<attrs.Size(); j++)
209  {
210  int attr = attrs[j];
211  MFEM_VERIFY(attr > 0, "Attribute number less than one!");
212  marker[attr-1] = 1;
213  }
214  }
215 }
216 
217 } // namespace common
218 
219 } // namespace mfem
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:917
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:849
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:908
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
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:291
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:153
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:942
int Dimension() const
Definition: mesh.hpp:911
int SpaceDimension() const
Definition: mesh.hpp:912
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:11057
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:674
int EulerNumber() const
Equals 1 + num_holes - num_loops.
Definition: mesh.hpp:905
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:843
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:946