GRINS-0.8.0
constrained_points.C
Go to the documentation of this file.
1 //-----------------------------------------------------------------------bl-
2 //--------------------------------------------------------------------------
3 //
4 // GRINS - General Reacting Incompressible Navier-Stokes
5 //
6 // Copyright (C) 2014-2017 Paul T. Bauman, Roy H. Stogner
7 // Copyright (C) 2010-2013 The PECOS Development Team
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the Version 2.1 GNU Lesser General
11 // Public License as published by the Free Software Foundation.
12 //
13 // This library is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // Lesser General Public License for more details.
17 //
18 // You should have received a copy of the GNU Lesser General Public
19 // License along with this library; if not, write to the Free Software
20 // Foundation, Inc. 51 Franklin Street, Fifth Floor,
21 // Boston, MA 02110-1301 USA
22 //
23 //-----------------------------------------------------------------------el-
24 
25 // This class
27 
28 // GRINS
29 #include "grins/builder_helper.h"
31 
32 // libMesh
33 #include "libmesh/dof_map.h"
34 #include "libmesh/getpot.h"
35 #include "libmesh/node.h"
36 #include "libmesh/system.h"
37 
38 namespace GRINS
39 {
41  libMesh::System& system )
42  : _sys(system)
43  {
44  std::vector<std::string> var_sections;
45 
46  BuilderHelper::parse_var_sections_vector(input, var_sections);
47 
48  for (std::vector<std::string>::const_iterator
49  it = var_sections.begin(); it != var_sections.end(); ++it)
50  {
51  const std::string & section_name = *it;
52  const std::string constraint_section =
54  section_name + '/' + "ConstrainedPoints";
55 
56  if (input.have_section (constraint_section))
57  {
58  std::vector<std::string> point_sections =
59  input.get_subsection_names(constraint_section);
60 
61  for (std::vector<std::string>::const_iterator
62  pt_it = point_sections.begin();
63  pt_it != point_sections.end(); ++pt_it)
64  {
66 
67  pt.name = *pt_it;
68 
69  const std::string point_section =
70  constraint_section + '/' + pt.name;
71 
72  // By default we constrain the first variable in our
73  // variable names list
74  std::string constrained_var_name =
76  section_name + '/' + "names", std::string(),
77  0);
78 
79  // But we can override that
80  constrained_var_name =
81  input(point_section + "/constraint_var",
82  constrained_var_name);
83 
84  pt.var = _sys.variable_number(constrained_var_name);
85 
86  const unsigned int constrained_loc_dim =
87  input.vector_variable_size(point_section+"/constraint_location");
88  libmesh_assert_less_equal(constrained_loc_dim, LIBMESH_DIM);
89  for (unsigned int d=0; d != constrained_loc_dim; ++d)
90  pt(d) = input(point_section+"/constraint_location", 0.0, d );
91 
92  const unsigned int n_constraining_points_x =
93  input.vector_variable_size(point_section+"/constraining_points_x");
94 
95 #ifndef NDEBUG
96  const unsigned int n_constraining_points_y =
97  input.vector_variable_size(point_section+"/constraining_points_y");
98  if (n_constraining_points_y)
99  libmesh_assert_equal_to(n_constraining_points_x, n_constraining_points_y);
100 
101  const unsigned int n_constraining_points_z =
102  input.vector_variable_size(point_section+"/constraining_points_z");
103  if (n_constraining_points_z)
104  libmesh_assert_equal_to(n_constraining_points_x, n_constraining_points_z);
105 
106  const unsigned int n_constraining_points_var =
107  input.vector_variable_size(point_section+"/constraining_points_var");
108  libmesh_assert_equal_to(n_constraining_points_x, n_constraining_points_var);
109 
110  const unsigned int n_constraining_points_coeff =
111  input.vector_variable_size(point_section+"/constraining_points_coeff");
112  if (n_constraining_points_coeff)
113  libmesh_assert_equal_to(n_constraining_points_x, n_constraining_points_coeff);
114 #endif // !NDEBUG
115 
116  for (unsigned int n=0; n != n_constraining_points_x; ++n)
117  {
118  ConstrainingPoint constraining_pt;
119  constraining_pt(0) = input(point_section+"/constraining_points_x", 0.0, n);
120  if (LIBMESH_DIM > 1)
121  constraining_pt(1) = input(point_section+"/constraining_points_y", 0.0, n);
122  if (LIBMESH_DIM > 2)
123  constraining_pt(2) = input(point_section+"/constraining_points_z", 0.0, n);
124 
125  constraining_pt.coeff =
126  input(point_section+"/constraining_points_coeff", libMesh::Number(0), n);
127 
128  const std::string var_name
129  (input(point_section+"/constraining_points_var", std::string(), n));
130 
131  constraining_pt.var = _sys.variable_number(var_name);
132 
133  pt.constrainers.push_back(constraining_pt);
134  }
135 
136  pt.rhs = input(point_section+"/constraint_rhs", 0.0);
137 
138  pt.forbid_overwrite = input(point_section+"/forbid_overwrite", true);
139 
140  _constrained_pts.push_back(pt);
141  }
142  }
143  }
144  }
145 
147  {
148  libMesh::MeshBase & mesh = _sys.get_mesh();
149  const unsigned int sys_num = _sys.number();
150 
151  libMesh::UniquePtr<libMesh::PointLocatorBase> locator =
152  mesh.sub_point_locator();
153 
154  for (std::vector<ConstrainedPoint>::const_iterator
155  it = _constrained_pts.begin(); it != _constrained_pts.end(); ++it)
156  {
157  const ConstrainedPoint & constrained_pt = *it;
158  const libMesh::Point & pt = constrained_pt;
159 
160  // Find the node at the constrained point
161  const libMesh::Node *constrained_node = locator->locate_node(pt);
162 
163  // If we're on a distributed mesh, not every processor might have found
164  // the node, but someone needs to have found it.
165  {
166  bool found_node = constrained_node;
167  _sys.comm().max(found_node);
168  if (!found_node)
169  libmesh_error_msg("Failed to find a Node at point " << pt);
170  }
171 
172  // We'll only need to add constraints on processors
173  // which have found the node. But those processors may not be
174  // able to find all the constraining nodes, so all processors
175  // need to stay in play here.
176 
177  const libMesh::dof_id_type constrained_dof =
178  constrained_node ? constrained_node->dof_number
179  (sys_num, constrained_pt.var, 0) : 0;
180 
181  libmesh_assert
182  (_sys.comm().semiverify (constrained_node ? &constrained_dof
183  : libmesh_nullptr));
184 
185  libMesh::DofConstraintRow constraint_row;
186 
187  for (std::vector<ConstrainingPoint>::const_iterator
188  jt = constrained_pt.constrainers.begin();
189  jt != constrained_pt.constrainers.end(); ++jt)
190  {
191  const ConstrainingPoint & constraining_pt = *jt;
192  const libMesh::Point & pt2 = constraining_pt;
193 
194  const libMesh::Node *constraining_node =
195  locator->locate_node(pt2);
196 
197  // If we're on a distributed mesh, not every processor
198  // might have found the node, but someone needs to have
199  // found it.
200  {
201  bool found_node = constraining_node;
202  _sys.comm().max(found_node);
203  if (!found_node)
204  libmesh_error_msg
205  ("Failed to find a Node at point " << pt);
206  }
207 
208  libMesh::dof_id_type constraining_dof =
209  constraining_node ? constraining_node->dof_number
210  (sys_num, constraining_pt.var, 0) : 0;
211 
212  libmesh_assert
213  (_sys.comm().semiverify (constraining_node ?
214  &constraining_dof :
215  libmesh_nullptr));
216 
217  _sys.comm().max(constraining_dof);
218 
219  if (constrained_node)
220  constraint_row[constraining_dof] = constraining_pt.coeff;
221  }
222 
223  if (constrained_node)
224  _sys.get_dof_map().add_constraint_row
225  (constrained_dof, constraint_row, constrained_pt.rhs,
226  constrained_pt.forbid_overwrite);
227  }
228  }
229 } // end namespace GRINS
std::vector< ConstrainingPoint > constrainers
static std::string variables_section()
Helper function to encapsualte the overall [Variables] section name.
ConstrainedPoints(const GetPot &input, libMesh::System &system)
static void parse_var_sections_vector(const GetPot &input, std::vector< std::string > &sections)
The same as parse_var_sections, except the result is returned in an std::vector.
GRINS namespace.
std::vector< ConstrainedPoint > _constrained_pts

Generated on Tue Dec 19 2017 12:47:27 for GRINS-0.8.0 by  doxygen 1.8.9.1