MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
mesh_extras.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 "mesh_extras.hpp"
13
14using namespace std;
15
16namespace mfem
17{
18
19namespace common
20{
21
22ElementMeshStream::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
139void
140MergeMeshNodes(Mesh * mesh, int logging)
141{
142 int dim = mesh->Dimension();
143 int sdim = mesh->SpaceDimension();
144
145 real_t 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 real_t 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
216void 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
237void 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 real_t x = pos(0), y = pos(1), z = dim == 3 ? pos(2) : 0;
244 real_t X, Y, Z;
245
246 X = x;
247
248 int layer = x*6.0;
249 real_t 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::Max
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
mfem::Array::SetSize
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
mfem::Array::Size
int Size() const
Return the logical size of the array.
Definition array.hpp:144
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:29
mfem::Element::GetVertices
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
mfem::Element::Type
Type
Constants for the classes derived from Element.
Definition element.hpp:41
mfem::Element::GetNVertices
virtual int GetNVertices() const =0
mfem::IntegrationPoint
Class for integration point with weight.
Definition intrules.hpp:35
mfem::Mesh
Mesh data type.
Definition mesh.hpp:56
mfem::Mesh::EulerNumber
int EulerNumber() const
Equals 1 + num_holes - num_loops.
Definition mesh.hpp:1166
mfem::Mesh::EulerNumber2D
int EulerNumber2D() const
Equals 1 - num_holes.
Definition mesh.hpp:1169
mfem::Mesh::GetElement
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
mfem::Mesh::Dimension
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
mfem::Mesh::GetBdrElement
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1298
mfem::Mesh::SpaceDimension
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
mfem::Mesh::GetCharacteristics
void GetCharacteristics(real_t &h_min, real_t &h_max, real_t &kappa_min, real_t &kappa_max, Vector *Vh=NULL, Vector *Vk=NULL)
Definition mesh.cpp:202
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:1223
mfem::Mesh::GetNBE
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1229
mfem::Mesh::GetVertex
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1265
mfem::Mesh::RemoveUnusedVertices
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition mesh.cpp:12872
mfem::Vector
Vector data type.
Definition vector.hpp:80
mfem::Vector::Norml2
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:538
mfem::common::KershawTransformation::right
real_t right(const real_t eps, const real_t x)
Definition mesh_extras.hpp:79
mfem::common::KershawTransformation::step
real_t step(const real_t a, const real_t b, real_t x)
Definition mesh_extras.hpp:92
mfem::common::KershawTransformation::Eval
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Definition mesh_extras.cpp:237
mfem::common::KershawTransformation::left
real_t left(const real_t eps, const real_t x)
Definition mesh_extras.hpp:85
dim
int dim
Definition ex24.cpp:53
mesh_extras.hpp
mfem::common::MergeMeshNodes
void MergeMeshNodes(Mesh *mesh, int logging)
Merges vertices which lie at the same location.
Definition mesh_extras.cpp:140
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::mfem_error
void mfem_error(const char *msg)
Definition error.cpp:154
mfem::add
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition vector.cpp:316
mfem::real_t
float real_t
Definition config.hpp:43
doxygen 1.11.0