libMesh
meshbcid.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // Open the mesh named in command line arguments,
19 // apply the named tests on every boundary side and normal vector,
20 // and give the specified boundary condition id to each side
21 // that passes the tests.
22 
23 #include <limits>
24 #include <string>
25 
26 #include "libmesh/libmesh.h"
27 
28 #include "libmesh/boundary_info.h"
29 #include "libmesh/elem.h"
30 #include "libmesh/fe.h"
31 #include "libmesh/getpot.h"
32 #include "libmesh/mesh.h"
33 #include "libmesh/point.h"
34 #include "libmesh/quadrature_gauss.h"
35 
36 using namespace libMesh;
37 
38 void usage_error(const char * progname)
39 {
40  libMesh::out << "Usage: " << progname
41  << " --dim d --input inputmesh --output outputmesh --newbcid idnum --tests --moretests"
42  << std::endl;
43 
44  exit(1);
45 }
46 
47 int main(int argc, char ** argv)
48 {
49  LibMeshInit init(argc, argv);
50 
51  GetPot cl(argc, argv);
52 
53  unsigned char dim = -1;
54  if (!cl.search("--dim"))
55  {
56  libMesh::err << "No --dim argument found!" << std::endl;
57  usage_error(argv[0]);
58  }
59  dim = cl.next(dim);
60 
61  Mesh mesh(init.comm(), dim);
62 
63  if (!cl.search("--input"))
64  {
65  libMesh::err << "No --input argument found!" << std::endl;
66  usage_error(argv[0]);
67  }
68  const char * meshname = cl.next("mesh.xda");
69 
70  mesh.read(meshname);
71  libMesh::out << "Loaded mesh " << meshname << std::endl;
72 
73  if (!cl.search("--newbcid"))
74  {
75  libMesh::err << "No --bcid argument found!" << std::endl;
76  usage_error(argv[0]);
77  }
78  boundary_id_type bcid = 0;
79  bcid = cl.next(bcid);
80 
93 
94  if (cl.search("--minnormalx"))
95  minnormal(0) = cl.next(minnormal(0));
96  if (cl.search("--minnormalx"))
97  minnormal(0) = cl.next(minnormal(0));
98  if (cl.search("--maxnormalx"))
99  maxnormal(0) = cl.next(maxnormal(0));
100  if (cl.search("--minnormaly"))
101  minnormal(1) = cl.next(minnormal(1));
102  if (cl.search("--maxnormaly"))
103  maxnormal(1) = cl.next(maxnormal(1));
104  if (cl.search("--minnormalz"))
105  minnormal(2) = cl.next(minnormal(2));
106  if (cl.search("--maxnormalz"))
107  maxnormal(2) = cl.next(maxnormal(2));
108 
109  if (cl.search("--minpointx"))
110  minpoint(0) = cl.next(minpoint(0));
111  if (cl.search("--maxpointx"))
112  maxpoint(0) = cl.next(maxpoint(0));
113  if (cl.search("--minpointy"))
114  minpoint(1) = cl.next(minpoint(1));
115  if (cl.search("--maxpointy"))
116  maxpoint(1) = cl.next(maxpoint(1));
117  if (cl.search("--minpointz"))
118  minpoint(2) = cl.next(minpoint(2));
119  if (cl.search("--maxpointz"))
120  maxpoint(2) = cl.next(maxpoint(2));
121 
122  libMesh::out << "min point = " << minpoint << std::endl;
123  libMesh::out << "max point = " << maxpoint << std::endl;
124  libMesh::out << "min normal = " << minnormal << std::endl;
125  libMesh::out << "max normal = " << maxnormal << std::endl;
126 
127  bool matcholdbcid = false;
128  boundary_id_type oldbcid = 0;
129  if (cl.search("--oldbcid"))
130  {
131  matcholdbcid = true;
132  oldbcid = cl.next(oldbcid);
133  if (oldbcid < 0)
134  oldbcid = BoundaryInfo::invalid_id;
135  }
136 
138  QGauss qface(dim-1, CONSTANT);
139  fe->attach_quadrature_rule(&qface);
140  const std::vector<Point> & face_points = fe->get_xyz();
141  const std::vector<Point> & face_normals = fe->get_normals();
142 
143  for (auto & elem : mesh.element_ptr_range())
144  {
145  unsigned int n_sides = elem->n_sides();
146 
147  // Container to catch ids handed back from BoundaryInfo
148  std::vector<boundary_id_type> ids;
149 
150  for (unsigned short s=0; s != n_sides; ++s)
151  {
152  if (elem->neighbor_ptr(s))
153  continue;
154 
155  fe->reinit(elem,s);
156  const Point & p = face_points[0];
157  const Point & n = face_normals[0];
158 
159  //libMesh::out << "elem = " << elem->id() << std::endl;
160  //libMesh::out << "centroid = " << elem->centroid() << std::endl;
161  //libMesh::out << "p = " << p << std::endl;
162  //libMesh::out << "n = " << n << std::endl;
163 
164  if (p(0) > minpoint(0) && p(0) < maxpoint(0) &&
165  p(1) > minpoint(1) && p(1) < maxpoint(1) &&
166  p(2) > minpoint(2) && p(2) < maxpoint(2) &&
167  n(0) > minnormal(0) && n(0) < maxnormal(0) &&
168  n(1) > minnormal(1) && n(1) < maxnormal(1) &&
169  n(2) > minnormal(2) && n(2) < maxnormal(2))
170  {
171  // Get the list of boundary ids for this side
172  mesh.get_boundary_info().boundary_ids(elem, s, ids);
173 
174  // There should be at most one value present, otherwise the
175  // logic here won't work.
176  libmesh_assert(ids.size() <= 1);
177 
178  // A convenient name for the side's ID.
179  boundary_id_type b_id = ids.empty() ? BoundaryInfo::invalid_id : ids[0];
180 
181  if (matcholdbcid && b_id != oldbcid)
182  continue;
183 
185  mesh.get_boundary_info().add_side(elem, s, bcid);
186  //libMesh::out << "Set element " << elem->id() << " side " << s <<
187  // " to boundary " << bcid << std::endl;
188  }
189  }
190  }
191 
192  // We might have removed *every* instance of a given id, and if that
193  // happened then we should make sure that file formats which write
194  // out id sets do not write out the removed id.
196 
197  std::string outputname;
198  if (cl.search("--output"))
199  {
200  outputname = cl.next("mesh.xda");
201  }
202  else
203  {
204  outputname = "new.";
205  outputname += meshname;
206  }
207 
208 
209  mesh.write(outputname.c_str());
210  libMesh::out << "Wrote mesh " << outputname << std::endl;
211 
212  return 0;
213 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
OStreamProxy err
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:117
unsigned int dim
MeshBase & mesh
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
std::vector< boundary_id_type > boundary_ids(const Node *node) const
int8_t boundary_id_type
Definition: id_types.h:51
static const boundary_id_type invalid_id
Number used for internal use.
virtual SimpleRange< element_iterator > element_ptr_range()=0
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
virtual void write(const std::string &name)=0
void regenerate_id_sets()
Clears and regenerates the cached sets of ids.
void remove_side(const Elem *elem, const unsigned short int side)
Removes all boundary conditions associated with side side of element elem, if any exist...
OStreamProxy out
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
Add side side of element number elem with boundary id id to the boundary information data structure...
This class implements specific orders of Gauss quadrature.
int main(int argc, char **argv)
Definition: meshbcid.C:47
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
virtual void read(const std::string &name, void *mesh_data=libmesh_nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
void usage_error(const char *progname)
Definition: meshbcid.C:38
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.