MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
trimmer.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// ------------------------------------------------------------------------
13// Trimmer Miniapp: Trim away elements according to their attribute numbers
14// ------------------------------------------------------------------------
15//
16// This miniapp creates a new mesh consisting of all the elements not possessing
17// a given set of attribute numbers. The new boundary elements are created with
18// boundary attribute numbers related to the trimmed elements' attribute
19// numbers.
20//
21// By default the new boundary elements will have new attribute numbers so as
22// not to interfere with existing boundaries. For example, consider a mesh with
23// attributes given by:
24//
25// attributes = {a1, a2, a3, a4, a5, a6, ..., amax}
26// bdr_attributes = {b1, b2, ..., bmax}
27//
28// If we trim away elements with attributes a2 and a4 the new mesh will have
29// attributes:
30//
31// attributes: {a1, a3, a5, a6, ..., amax}
32// bdr_attributes = {b1, b2, ..., bmax, bmax + a2, bmax + a4}
33//
34// The user has the option of providing new attribute numbers for each group of
35// elements to be trimmed. In this case the new boundary elements may have the
36// same attribute numbers as existing boundary elements.
37//
38// The resulting mesh is displayed with GLVis (unless explicitly disabled) and
39// is also written to the file "trimmer.mesh"
40//
41// Compile with: make trimmer
42//
43// Sample runs: trimmer -a '2' -b '2'
44// trimmer -m ../../data/beam-hex.mesh -a '2'
45// trimmer -m ../../data/beam-hex.mesh -a '2' -b '2'
46
47#include "mfem.hpp"
48#include <fstream>
49#include <iostream>
50
51using namespace std;
52using namespace mfem;
53
54int main(int argc, char *argv[])
55{
56 // Parse command-line options.
57 const char *mesh_file = "../../data/beam-tet.vtk";
58 Array<int> attr;
59 Array<int> bdr_attr;
60 bool visualization = 1;
61
62 OptionsParser args(argc, argv);
63 args.AddOption(&mesh_file, "-m", "--mesh",
64 "Mesh file to use.");
65 args.AddOption(&attr, "-a", "--attr",
66 "Set of attributes to remove from the mesh.");
67 args.AddOption(&bdr_attr, "-b", "--bdr-attr",
68 "Set of attributes to assign to the new boundary elements.");
69 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
70 "--no-visualization",
71 "Enable or disable GLVis visualization.");
72 args.Parse();
73 if (!args.Good())
74 {
75 args.PrintUsage(cout);
76 return 1;
77 }
78 args.PrintOptions(cout);
79
80 Mesh mesh(mesh_file, 0, 0);
81
82 int max_attr = mesh.attributes.Max();
83 int max_bdr_attr = mesh.bdr_attributes.Max();
84
85 if (bdr_attr.Size() == 0)
86 {
87 bdr_attr.SetSize(attr.Size());
88 for (int i=0; i<attr.Size(); i++)
89 {
90 bdr_attr[i] = max_bdr_attr + attr[i];
91 }
92 }
93 MFEM_VERIFY(attr.Size() == bdr_attr.Size(),
94 "Size mismatch in attribute arguments.");
95
96 Array<int> marker(max_attr);
97 Array<int> attr_inv(max_attr);
98 marker = 0;
99 attr_inv = 0;
100 for (int i=0; i<attr.Size(); i++)
101 {
102 marker[attr[i]-1] = 1;
103 attr_inv[attr[i]-1] = i;
104 }
105
106 // Count the number of elements in the final mesh
107 int num_elements = 0;
108 for (int e=0; e<mesh.GetNE(); e++)
109 {
110 int elem_attr = mesh.GetElement(e)->GetAttribute();
111 if (!marker[elem_attr-1]) { num_elements++; }
112 }
113
114 // Count the number of boundary elements in the final mesh
115 int num_bdr_elements = 0;
116 for (int f=0; f<mesh.GetNumFaces(); f++)
117 {
118 int e1 = -1, e2 = -1;
119 mesh.GetFaceElements(f, &e1, &e2);
120
121 int a1 = 0, a2 = 0;
122 if (e1 >= 0) { a1 = mesh.GetElement(e1)->GetAttribute(); }
123 if (e2 >= 0) { a2 = mesh.GetElement(e2)->GetAttribute(); }
124
125 if (a1 == 0 || a2 == 0)
126 {
127 if (a1 == 0 && !marker[a2-1]) { num_bdr_elements++; }
128 else if (a2 == 0 && !marker[a1-1]) { num_bdr_elements++; }
129 }
130 else
131 {
132 if (marker[a1-1] && !marker[a2-1]) { num_bdr_elements++; }
133 else if (!marker[a1-1] && marker[a2-1]) { num_bdr_elements++; }
134 }
135 }
136
137 cout << "Number of Elements: " << mesh.GetNE() << " -> "
138 << num_elements << endl;
139 cout << "Number of Boundary Elements: " << mesh.GetNBE() << " -> "
140 << num_bdr_elements << endl;
141
142 Mesh trimmed_mesh(mesh.Dimension(), mesh.GetNV(),
143 num_elements, num_bdr_elements, mesh.SpaceDimension());
144
145 // Copy vertices
146 for (int v=0; v<mesh.GetNV(); v++)
147 {
148 trimmed_mesh.AddVertex(mesh.GetVertex(v));
149 }
150
151 // Copy elements
152 for (int e=0; e<mesh.GetNE(); e++)
153 {
154 Element * el = mesh.GetElement(e);
155 int elem_attr = el->GetAttribute();
156 if (!marker[elem_attr-1])
157 {
158 Element * nel = mesh.NewElement(el->GetGeometryType());
159 nel->SetAttribute(elem_attr);
160 nel->SetVertices(el->GetVertices());
161 trimmed_mesh.AddElement(nel);
162 }
163 }
164
165 // Copy selected boundary elements
166 for (int be=0; be<mesh.GetNBE(); be++)
167 {
168 int e, info;
169 mesh.GetBdrElementAdjacentElement(be, e, info);
170
171 int elem_attr = mesh.GetElement(e)->GetAttribute();
172 if (!marker[elem_attr-1])
173 {
174 Element * nbel = mesh.GetBdrElement(be)->Duplicate(&trimmed_mesh);
175 trimmed_mesh.AddBdrElement(nbel);
176 }
177 }
178
179 // Create new boundary elements
180 for (int f=0; f<mesh.GetNumFaces(); f++)
181 {
182 int e1 = -1, e2 = -1;
183 mesh.GetFaceElements(f, &e1, &e2);
184
185 int i1 = -1, i2 = -1;
186 mesh.GetFaceInfos(f, &i1, &i2);
187
188 int a1 = 0, a2 = 0;
189 if (e1 >= 0) { a1 = mesh.GetElement(e1)->GetAttribute(); }
190 if (e2 >= 0) { a2 = mesh.GetElement(e2)->GetAttribute(); }
191
192 if (a1 != 0 && a2 != 0)
193 {
194 if (marker[a1-1] && !marker[a2-1])
195 {
196 Element * bel = (mesh.Dimension() == 1) ?
197 (Element*)new Point(&f) :
198 mesh.GetFace(f)->Duplicate(&trimmed_mesh);
199 bel->SetAttribute(bdr_attr[attr_inv[a1-1]]);
200 trimmed_mesh.AddBdrElement(bel);
201 }
202 else if (!marker[a1-1] && marker[a2-1])
203 {
204 Element * bel = (mesh.Dimension() == 1) ?
205 (Element*)new Point(&f) :
206 mesh.GetFace(f)->Duplicate(&trimmed_mesh);
207 bel->SetAttribute(bdr_attr[attr_inv[a2-1]]);
208 trimmed_mesh.AddBdrElement(bel);
209 }
210 }
211 }
212
213 trimmed_mesh.FinalizeTopology();
214 trimmed_mesh.Finalize();
215 trimmed_mesh.RemoveUnusedVertices();
216
217 // Check for curved or discontinuous mesh
218 if (mesh.GetNodes())
219 {
220 // Extract Nodes GridFunction and determine its type
221 const GridFunction * Nodes = mesh.GetNodes();
222 const FiniteElementSpace * fes = Nodes->FESpace();
223
224 Ordering::Type ordering = fes->GetOrdering();
225 int order = fes->FEColl()->GetOrder();
226 int sdim = mesh.SpaceDimension();
227 bool discont =
228 dynamic_cast<const L2_FECollection*>(fes->FEColl()) != NULL;
229
230 // Set curvature of the same type as original mesh
231 trimmed_mesh.SetCurvature(order, discont, sdim, ordering);
232
233 const FiniteElementSpace * trimmed_fes = trimmed_mesh.GetNodalFESpace();
234 GridFunction * trimmed_nodes = trimmed_mesh.GetNodes();
235
236 Array<int> vdofs;
237 Array<int> trimmed_vdofs;
238 Vector loc_vec;
239
240 // Copy nodes to trimmed mesh
241 int te = 0;
242 for (int e = 0; e < mesh.GetNE(); e++)
243 {
244 Element * el = mesh.GetElement(e);
245 int elem_attr = el->GetAttribute();
246 if (!marker[elem_attr-1])
247 {
248 fes->GetElementVDofs(e, vdofs);
249 Nodes->GetSubVector(vdofs, loc_vec);
250
251 trimmed_fes->GetElementVDofs(te, trimmed_vdofs);
252 trimmed_nodes->SetSubVector(trimmed_vdofs, loc_vec);
253 te++;
254 }
255 }
256 }
257
258 // Save the final mesh
259 ofstream mesh_ofs("trimmer.mesh");
260 mesh_ofs.precision(8);
261 trimmed_mesh.Print(mesh_ofs);
262
263 if (visualization)
264 {
265 // GLVis server to visualize to
266 char vishost[] = "localhost";
267 int visport = 19916;
268 socketstream sol_sock(vishost, visport);
269 sol_sock.precision(8);
270 sol_sock << "mesh\n" << trimmed_mesh << flush;
271 }
272}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:697
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
Abstract data type element.
Definition: element.hpp:29
Geometry::Type GetGeometryType() const
Definition: element.hpp:52
virtual Element * Duplicate(Mesh *m) const =0
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
void SetAttribute(const int attr)
Set element's attribute.
Definition: element.hpp:58
int GetAttribute() const
Return element's attribute.
Definition: element.hpp:55
virtual void SetVertices(const Array< int > &v)=0
Set the indices defining the vertices.
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:225
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i'th element. The returned indices are offsets into an ...
Definition: fespace.cpp:287
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:725
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:727
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:696
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
Mesh data type.
Definition: mesh.hpp:56
Element * NewElement(int geom)
Definition: mesh.cpp:4401
int AddBdrElement(Element *elem)
Definition: mesh.cpp:2028
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:282
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition: mesh.cpp:1439
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:6250
const FiniteElementSpace * GetNodalFESpace() const
Definition: mesh.cpp:6206
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition: mesh.hpp:1283
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:3135
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition: mesh.hpp:2288
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition: mesh.cpp:1658
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
const Element * GetFace(int i) const
Return pointer to the i'th face element object.
Definition: mesh.hpp:1310
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition: mesh.hpp:1298
int AddElement(Element *elem)
Definition: mesh.cpp:2021
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1433
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1163
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8973
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:1223
void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
For the given boundary element, bdr_el, return its adjacent element and its info, i....
Definition: mesh.cpp:7271
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1229
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition: mesh.cpp:3241
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:6211
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:280
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition: mesh.hpp:1265
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:12872
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Type
Ordering methods:
Definition: fespace.hpp:34
Data type point element.
Definition: point.hpp:23
Vector data type.
Definition: vector.hpp:80
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:604
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:578
int visport
char vishost[]
int main()
real_t f(const Vector &p)