MFEM  v4.6.0
Finite element discretization library
trimmer.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 
51 using namespace std;
52 using namespace mfem;
53 
54 int 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 }
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
int visport
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
virtual Element * Duplicate(Mesh *m) const =0
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5668
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:225
int GetAttribute() const
Return element&#39;s attribute.
Definition: element.hpp:55
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
virtual void SetVertices(const int *ind)
Set the indices the element according to the input.
Definition: element.cpp:17
const Element * GetElement(int i) const
Return pointer to the i&#39;th element object.
Definition: mesh.hpp:1143
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
STL namespace.
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:579
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1089
Element * NewElement(int geom)
Definition: mesh.cpp:3901
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1402
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:1083
int AddVertex(double x, double y=0.0, double z=0.0)
Definition: mesh.cpp:1626
Geometry::Type GetGeometryType() const
Definition: element.hpp:52
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
Type
Ordering methods:
Definition: fespace.hpp:33
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1125
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition: mesh.cpp:1408
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
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 &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
int main(int argc, char *argv[])
Definition: trimmer.cpp:54
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:721
const Element * GetFace(int i) const
Definition: mesh.hpp:1167
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:58
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Vector data type.
Definition: vector.hpp:58
Data type point element.
Definition: point.hpp:22
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8302
const Element * GetBdrElement(int i) const
Return pointer to the i&#39;th boundary element object.
Definition: mesh.hpp:1158
Abstract data type element.
Definition: element.hpp:28
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:273
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:6600
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:327
double f(const Vector &p)