MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
mesh_extras.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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 case Element::PYRAMID:
133 *this << "dimension" << endl << 3 << endl
134 << "elements" << endl << 1 << endl
135 << "1 7 0 1 2 3 4" << endl
136 << "boundary" << endl << 5 << endl
137 << "1 3 3 2 1 0" << endl
138 << "1 2 0 1 4" << endl
139 << "1 2 1 2 4" << endl
140 << "1 2 3 4 2" << endl
141 << "1 2 0 4 3" << endl
142 << "vertices" << endl
143 << "5" << endl
144 << "3" << endl
145 << "0 0 0" << endl
146 << "1 0 0" << endl
147 << "1 1 0" << endl
148 << "0 1 0" << endl
149 << "0 0 1" << endl;
150 break;
151 default:
152 mfem_error("Invalid element type!");
153 break;
154 }
155
156}
157
158void
159MergeMeshNodes(Mesh * mesh, int logging)
160{
161 int dim = mesh->Dimension();
162 int sdim = mesh->SpaceDimension();
163
164 real_t h_min, h_max, k_min, k_max;
165 mesh->GetCharacteristics(h_min, h_max, k_min, k_max);
166
167 // Set tolerance for merging vertices
168 real_t tol = 1.0e-8 * h_min;
169
170 if ( logging > 0 )
171 cout << "Euler Number of Initial Mesh: "
172 << ((dim==3)?mesh->EulerNumber() :
173 ((dim==2)?mesh->EulerNumber2D() :
174 mesh->GetNV() - mesh->GetNE())) << endl;
175
176 vector<int> v2v(mesh->GetNV());
177
178 Vector vd(sdim);
179
180 for (int i = 0; i < mesh->GetNV(); i++)
181 {
182 Vector vi(mesh->GetVertex(i), sdim);
183
184 v2v[i] = -1;
185
186 for (int j = 0; j < i; j++)
187 {
188 Vector vj(mesh->GetVertex(j), sdim);
189 add(vi, -1.0, vj, vd);
190
191 if ( vd.Norml2() < tol )
192 {
193 v2v[i] = j;
194 break;
195 }
196 }
197 if ( v2v[i] < 0 ) { v2v[i] = i; }
198 }
199
200 // renumber elements
201 for (int i = 0; i < mesh->GetNE(); i++)
202 {
203 Element *el = mesh->GetElement(i);
204 int *v = el->GetVertices();
205 int nv = el->GetNVertices();
206 for (int j = 0; j < nv; j++)
207 {
208 v[j] = v2v[v[j]];
209 }
210 }
211 // renumber boundary elements
212 for (int i = 0; i < mesh->GetNBE(); i++)
213 {
214 Element *el = mesh->GetBdrElement(i);
215 int *v = el->GetVertices();
216 int nv = el->GetNVertices();
217 for (int j = 0; j < nv; j++)
218 {
219 v[j] = v2v[v[j]];
220 }
221 }
222
223 mesh->RemoveUnusedVertices();
224
225 if ( logging > 0 )
226 {
227 cout << "Euler Number of Final Mesh: "
228 << ((dim==3) ? mesh->EulerNumber() :
229 ((dim==2) ? mesh->EulerNumber2D() :
230 mesh->GetNV() - mesh->GetNE()))
231 << endl;
232 }
233}
234
235void AttrToMarker(int max_attr, const Array<int> &attrs, Array<int> &marker)
236{
237 MFEM_ASSERT(attrs.Max() <= max_attr, "Invalid attribute number present.");
238
239 marker.SetSize(max_attr);
240 if (attrs.Size() == 1 && attrs[0] == -1)
241 {
242 marker = 1;
243 }
244 else
245 {
246 marker = 0;
247 for (int j=0; j<attrs.Size(); j++)
248 {
249 int attr = attrs[j];
250 MFEM_VERIFY(attr > 0, "Attribute number less than one!");
251 marker[attr-1] = 1;
252 }
253 }
254}
255
256void KershawTransformation::Eval(Vector &V, ElementTransformation &T,
257 const IntegrationPoint &ip)
258{
259 V = 0.0;
260 Vector pos(dim);
261 T.Transform(ip, pos);
262 real_t x = pos(0), y = pos(1), z = dim == 3 ? pos(2) : 0;
263 real_t X, Y, Z;
264
265 X = x;
266
267 int layer = x*6.0;
268 real_t lambda = (x-layer/6.0)*6;
269
270 // The x-range is split in 6 layers going from left-to-left, left-to-right,
271 // right-to-left (2 layers), left-to-right and right-to-right yz-faces.
272 switch (layer)
273 {
274 case 0:
275 Y = left(epsy, y);
276 Z = left(epsz, z);
277 break;
278 case 1:
279 case 4:
280 Y = step(left(epsy, y), right(epsy, y), lambda);
281 Z = step(left(epsz, z), right(epsz, z), lambda);
282 break;
283 case 2:
284 Y = step(right(epsy, y), left(epsy, y), lambda/2);
285 Z = step(right(epsz, z), left(epsz, z), lambda/2);
286 break;
287 case 3:
288 Y = step(right(epsy, y), left(epsy, y), (1+lambda)/2);
289 Z = step(right(epsz, z), left(epsz, z), (1+lambda)/2);
290 break;
291 default:
292 Y = right(epsy, y);
293 Z = right(epsz, z);
294 break;
295 }
296 V.SetSize(dim);
297 V(0) = X;
298 V(1) = Y;
299 if (dim == 3) { V(2) = Z; }
300}
301
302void SpiralTransformation::Eval(Vector &V, ElementTransformation &T,
303 const IntegrationPoint &ip)
304{
305 Vector pos(dim);
306 T.Transform(ip, pos);
307 real_t x = pos(0), y = pos(1), z = dim == 3 ? pos(2) : 0;
308
309 real_t theta = 2.0*M_PI*turns*x;
310 real_t r_min = (0.5-0.5*width) + (gap+width)*turns*x;
311 real_t r_xyz = r_min + (width)*y;
312
313 V.SetSize(dim);
314 V(0) = r_xyz*std::cos(theta);
315 V(1) = r_xyz*std::sin(theta);
316 if (dim == 3) { V(2) = z*width + x*height; }
317}
318
319} // namespace common
320
321} // 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:758
mfem::Array::Size
int Size() const
Return the logical size of the array.
Definition array.hpp:147
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:64
mfem::Mesh::EulerNumber
int EulerNumber() const
Equals 1 + num_holes - num_loops.
Definition mesh.hpp:1222
mfem::Mesh::EulerNumber2D
int EulerNumber2D() const
Equals 1 - num_holes.
Definition mesh.hpp:1225
mfem::Mesh::GetElement
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1339
mfem::Mesh::GetNE
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1282
mfem::Mesh::Dimension
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1216
mfem::Mesh::GetBdrElement
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition mesh.hpp:1354
mfem::Mesh::SpaceDimension
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1219
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:1279
mfem::Mesh::GetNBE
int GetNBE() const
Returns number of boundary elements.
Definition mesh.hpp:1285
mfem::Mesh::GetVertex
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1321
mfem::Mesh::RemoveUnusedVertices
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition mesh.cpp:13197
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::Vector::Norml2
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:931
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
mfem::common::KershawTransformation::right
real_t right(const real_t eps, const real_t x)
Definition mesh_extras.hpp:79
mfem::common::KershawTransformation::Eval
void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip) override
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Definition mesh_extras.cpp:256
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::left
real_t left(const real_t eps, const real_t x)
Definition mesh_extras.hpp:85
mfem::common::SpiralTransformation::Eval
void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip) override
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Definition mesh_extras.cpp:302
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:159
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:235
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:391
mfem::real_t
float real_t
Definition config.hpp:43
Generated on Wed Apr 9 2025 16:53:05 for MFEM by doxygen 1.11.0