MAST3
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
bracket3d.hpp
Go to the documentation of this file.
1 /*
2  * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit
3  * Copyright (C) 2013-2020 Manav Bhatia and MAST authors
4  *
5  * This library is free software; you can redistribute it and/or
6  * modify it under the terms of the GNU Lesser General Public
7  * License as published by the Free Software Foundation; either
8  * version 2.1 of the License, or (at your option) any later version.
9  *
10  * This library is distributed in the hope that it will be useful,
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13  * Lesser General Public License for more details.
14  *
15  * You should have received a copy of the GNU Lesser General Public
16  * License along with this library; if not, write to the Free Software
17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18  */
19 
20 #ifndef __mast_mesh_generation_bracket_3d_h__
21 #define __mast_mesh_generation_bracket_3d_h__
22 
23 // MAST includes
24 #include <mast/base/mast_data_types.h>
25 #include <mast/base/parameter_data.hpp>
26 #include <mast/util/getpot_wrapper.hpp>
27 #include <mast/optimization/design_parameter_vector.hpp>
28 
29 // libMesh includes
30 #include <libmesh/system.h>
31 #include <libmesh/unstructured_mesh.h>
32 #include <libmesh/fe_type.h>
33 #include <libmesh/string_to_enum.h>
34 #include <libmesh/mesh_generation.h>
35 #include <libmesh/elem.h>
36 #include <libmesh/node.h>
37 #include <libmesh/boundary_info.h>
38 #include <libmesh/dirichlet_boundaries.h>
39 #include <libmesh/zero_function.h>
40 #include <libmesh/mesh_refinement.h>
41 #include <libmesh/mesh_modification.h>
42 
43 
44 namespace MAST {
45 namespace Mesh {
46 namespace Generation {
47 
48 
49 struct Bracket3D {
50 
51  template <typename ScalarType>
52  class Pressure {
53  public:
54  Pressure(real_t p,
55  real_t l1,
56  real_t frac):
57  _p(p), _l1(l1), _frac(frac)
58  {}
59  virtual ~Pressure() {}
60 
61  template <typename ContextType>
62  inline ScalarType value(ContextType& c) const {
63  ScalarType v=(c.qp_location(0)>=_l1*(1.-_frac))?_p:0.;
64  return v;
65  }
66 
67  template <typename ContextType, typename ScalarFieldType>
68  inline ScalarType derivative(ContextType& c,
69  const ScalarFieldType& f) const {
70  return 0.;
71  }
72 
73  private:
74  real_t _p, _l1, _frac;
75  };
76 
77  static const uint_t dim = 3;
78  template <typename ScalarType>
79  using pressure_t = MAST::Mesh::Generation::Bracket3D::Pressure<ScalarType>;
80 
81  template <typename Context>
82  inline real_t
83  reference_volume(Context& c) {
84 
85  real_t
86  length = c.input("length", "length of domain along x-axis", 0.3),
87  height = c.input("height", "length of domain along y-axis", 0.3),
88  width = c.input("width", "length of domain along z-axis", 0.06);
89 
90  return length * height * width;
91  }
92 
93 
94  inline uint_t idx(const libMesh::ElemType type,
95  const uint_t nx,
96  const uint_t ny,
97  const uint_t i,
98  const uint_t j,
99  const uint_t k)
100  {
101  switch(type)
102  {
103  case libMesh::HEX8:
104  {
105  return i + (nx+1)*(j + k*(ny+1));
106  }
107 
108  case libMesh::HEX27:
109  {
110  return i + (2*nx+1)*(j + k*(2*ny+1));
111  }
112 
113  default:
114  Error(false, "Invalid element type");
115  }
116 
117  return libMesh::invalid_uint;
118  }
119 
120 
121  inline void build_cube(libMesh::UnstructuredMesh & mesh,
122  const uint_t nx,
123  const uint_t ny,
124  const uint_t nz,
125  const real_t length,
126  const real_t height,
127  const real_t width,
128  const libMesh::ElemType type) {
129 
130  Assert0(type == libMesh::HEX8 || type == libMesh::HEX27,
131  "Method only implemented for Hex8/Hex27");
132 
133  // Clear the mesh and start from scratch
134  mesh.clear();
135 
136  libMesh::BoundaryInfo & boundary_info = mesh.get_boundary_info();
137 
138  mesh.set_mesh_dimension(3);
139  mesh.set_spatial_dimension(3);
140  mesh.reserve_elem(nx*ny*nz);
141 
142  if (type == libMesh::HEX8)
143  mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
144  else if (type == libMesh::HEX27)
145  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
146 
147  real_t
148  xmax = length,
149  ymax = height,
150  zmax = width;
151 
152 
153 
154  std::map<uint_t, libMesh::Node*> nodes;
155 
156  // Build the nodes.
157  uint_t
158  node_id = 0,
159  n = 0;
160  switch (type)
161  {
162  case libMesh::HEX8: {
163 
164  for (uint_t k=0; k<=nz; k++)
165  for (uint_t j=0; j<=ny; j++)
166  for (uint_t i=0; i<=nx; i++) {
167  if ( i<=nx/10*4 || j>=ny/10*6) {
168  nodes[node_id] =
169  mesh.add_point(libMesh::Point(static_cast<real_t>(i)/static_cast<real_t>(nx)*length,
170  static_cast<real_t>(j)/static_cast<real_t>(ny)*height,
171  static_cast<real_t>(k)/static_cast<real_t>(nz)*width),
172  n++);
173  }
174  node_id++;
175  }
176 
177 
178  break;
179  }
180 
181  case libMesh::HEX27: {
182 
183  for (uint_t k=0; k<=(2*nz); k++)
184  for (uint_t j=0; j<=(2*ny); j++)
185  for (uint_t i=0; i<=(2*nx); i++) {
186  if ( i<=2*nx/10*4 || j>=2*ny/10*6) {
187  nodes[node_id] =
188  mesh.add_point(libMesh::Point(static_cast<real_t>(i)/static_cast<real_t>(2*nx)*length,
189  static_cast<real_t>(j)/static_cast<real_t>(2*ny)*height,
190  static_cast<real_t>(k)/static_cast<real_t>(2*nz)*width),
191  n++);
192  }
193  node_id++;
194  }
195 
196  break;
197  }
198 
199 
200  default:
201  Assert0(false, "ERROR: Unrecognized 3D element type.");
202  }
203 
204  // Build the elements.
205  uint_t
206  elem_id = 0;
207  switch (type)
208  {
209  case libMesh::HEX8:
210  {
211  for (uint_t k=0; k<nz; k++)
212  for (uint_t j=0; j<ny; j++)
213  for (uint_t i=0; i<nx; i++) {
214  if (i < nx*4/10 || j>=ny*6/10) {
215 
216  libMesh::Elem
217  *elem = libMesh::Elem::build(libMesh::HEX8).release();
218  elem->set_id(elem_id++);
219  mesh.add_elem(elem);
220 
221  elem->set_node(0) = nodes[idx(type,nx,ny,i,j,k) ];
222  elem->set_node(1) = nodes[idx(type,nx,ny,i+1,j,k) ];
223  elem->set_node(2) = nodes[idx(type,nx,ny,i+1,j+1,k) ];
224  elem->set_node(3) = nodes[idx(type,nx,ny,i,j+1,k) ];
225  elem->set_node(4) = nodes[idx(type,nx,ny,i,j,k+1) ];
226  elem->set_node(5) = nodes[idx(type,nx,ny,i+1,j,k+1) ];
227  elem->set_node(6) = nodes[idx(type,nx,ny,i+1,j+1,k+1)];
228  elem->set_node(7) = nodes[idx(type,nx,ny,i,j+1,k+1) ];
229 
230  if (k == 0)
231  boundary_info.add_side(elem, 0, 0);
232 
233  if (k == (nz-1))
234  boundary_info.add_side(elem, 5, 5);
235 
236  if (j == 0)
237  boundary_info.add_side(elem, 1, 1);
238 
239  if (j == (ny-1))
240  boundary_info.add_side(elem, 3, 3);
241 
242  if (i == 0)
243  boundary_info.add_side(elem, 4, 4);
244 
245  if (i == (nx-1))
246  boundary_info.add_side(elem, 2, 2);
247 
248  if (i==nx*3/10 && j<ny*6/10)
249  boundary_info.add_side(elem, 2, 6);
250 
251  if (j==ny*6/10 && i>nx*3/10)
252  boundary_info.add_side(elem, 1, 7);
253  }
254  }
255  break;
256  }
257 
258 
259  case libMesh::HEX27: {
260 
261  for (uint_t k=0; k<(2*nz); k += 2)
262  for (uint_t j=0; j<(2*ny); j += 2)
263  for (uint_t i=0; i<(2*nx); i += 2) {
264  if (i < 2*nx*4/10 || j>=2*ny*6/10) {
265 
266  libMesh::Elem
267  *elem = libMesh::Elem::build(libMesh::HEX27).release();
268  elem->set_id(elem_id++);
269  mesh.add_elem(elem);
270 
271  elem->set_node(0) = nodes[idx(type,nx,ny,i, j, k) ];
272  elem->set_node(1) = nodes[idx(type,nx,ny,i+2,j, k) ];
273  elem->set_node(2) = nodes[idx(type,nx,ny,i+2,j+2,k) ];
274  elem->set_node(3) = nodes[idx(type,nx,ny,i, j+2,k) ];
275  elem->set_node(4) = nodes[idx(type,nx,ny,i, j, k+2)];
276  elem->set_node(5) = nodes[idx(type,nx,ny,i+2,j, k+2)];
277  elem->set_node(6) = nodes[idx(type,nx,ny,i+2,j+2,k+2)];
278  elem->set_node(7) = nodes[idx(type,nx,ny,i, j+2,k+2)];
279  elem->set_node(8) = nodes[idx(type,nx,ny,i+1,j, k) ];
280  elem->set_node(9) = nodes[idx(type,nx,ny,i+2,j+1,k) ];
281  elem->set_node(10) = nodes[idx(type,nx,ny,i+1,j+2,k) ];
282  elem->set_node(11) = nodes[idx(type,nx,ny,i, j+1,k) ];
283  elem->set_node(12) = nodes[idx(type,nx,ny,i, j, k+1)];
284  elem->set_node(13) = nodes[idx(type,nx,ny,i+2,j, k+1)];
285  elem->set_node(14) = nodes[idx(type,nx,ny,i+2,j+2,k+1)];
286  elem->set_node(15) = nodes[idx(type,nx,ny,i, j+2,k+1)];
287  elem->set_node(16) = nodes[idx(type,nx,ny,i+1,j, k+2)];
288  elem->set_node(17) = nodes[idx(type,nx,ny,i+2,j+1,k+2)];
289  elem->set_node(18) = nodes[idx(type,nx,ny,i+1,j+2,k+2)];
290  elem->set_node(19) = nodes[idx(type,nx,ny,i, j+1,k+2)];
291 
292  elem->set_node(20) = nodes[idx(type,nx,ny,i+1,j+1,k) ];
293  elem->set_node(21) = nodes[idx(type,nx,ny,i+1,j, k+1)];
294  elem->set_node(22) = nodes[idx(type,nx,ny,i+2,j+1,k+1)];
295  elem->set_node(23) = nodes[idx(type,nx,ny,i+1,j+2,k+1)];
296  elem->set_node(24) = nodes[idx(type,nx,ny,i, j+1,k+1)];
297  elem->set_node(25) = nodes[idx(type,nx,ny,i+1,j+1,k+2)];
298  elem->set_node(26) = nodes[idx(type,nx,ny,i+1,j+1,k+1)];
299 
300  if (k == 0)
301  boundary_info.add_side(elem, 0, 0);
302 
303  if (k == 2*(nz-1))
304  boundary_info.add_side(elem, 5, 5);
305 
306  if (j == 0)
307  boundary_info.add_side(elem, 1, 1);
308 
309  if (j == 2*(ny-1))
310  boundary_info.add_side(elem, 3, 3);
311 
312  if (i == 0)
313  boundary_info.add_side(elem, 4, 4);
314 
315  if (i == 2*(nx-1))
316  boundary_info.add_side(elem, 2, 2);
317 
318  if (i==2*nx*3/10 && j<2*ny*6/10)
319  boundary_info.add_side(elem, 2, 6);
320 
321  if (j==2*ny*6/10 && i>2*nx*3/10)
322  boundary_info.add_side(elem, 1, 7);
323  }
324  }
325  break;
326  }
327 
328  default:
329  Assert0(false, "ERROR: Unrecognized 3D element type.");
330  }
331 
332  // Add sideset names to boundary info (Z axis out of the screen)
333  boundary_info.sideset_name(0) = "back";
334  boundary_info.sideset_name(1) = "bottom";
335  boundary_info.sideset_name(2) = "right";
336  boundary_info.sideset_name(3) = "top";
337  boundary_info.sideset_name(4) = "left";
338  boundary_info.sideset_name(5) = "front";
339  boundary_info.sideset_name(6) = "facing_right";
340  boundary_info.sideset_name(7) = "facing_down";
341 
342  // Add nodeset names to boundary info
343  boundary_info.nodeset_name(0) = "back";
344  boundary_info.nodeset_name(1) = "bottom";
345  boundary_info.nodeset_name(2) = "right";
346  boundary_info.nodeset_name(3) = "top";
347  boundary_info.nodeset_name(4) = "left";
348  boundary_info.nodeset_name(5) = "front";
349 
350  // Done building the mesh. Now prepare it for use.
351  mesh.prepare_for_use ();
352  }
353 
354 
355 
356  template <typename Context>
357  inline void
358  init_analysis_mesh(Context& c,
359  libMesh::UnstructuredMesh& mesh) {
360 
361  real_t
362  length = c.input("length", "length of domain along x-axis", 0.3),
363  height = c.input("height", "length of domain along y-axis", 0.3),
364  width = c.input("width", "length of domain along z-axis", 0.06);
365 
366  uint_t
367  nx_divs = c.input("nx_divs", "number of elements along x-axis", 10),
368  ny_divs = c.input("ny_divs", "number of elements along y-axis", 10),
369  nz_divs = c.input("nz_divs", "number of elements along z-axis", 2),
370  n_refine= c.input("n_uniform_refinement", "number of times the mesh is uniformly refined", 0);
371 
372  if (nx_divs%10 != 0 || ny_divs%10 != 0)
373  Error(false, "number of divisions in x and y must be multiples of 10");
374 
375  std::string
376  t = c.input("elem_type", "type of geometric element in the mesh", "hex8");
377 
378  libMesh::ElemType
379  e_type = libMesh::Utility::string_to_enum<libMesh::ElemType>(t);
380 
381  //
382  // if high order FE is used, libMesh requires atleast a second order
383  // geometric element.
384  //
385  if (c.sol_fe_order > 1 && e_type == libMesh::HEX8)
386  e_type = libMesh::HEX27;
387  else if (c.sol_fe_order > 1 && e_type == libMesh::TET4)
388  e_type = libMesh::TET10;
389 
390  //
391  // initialize the mesh with one element
392  //
393  build_cube(mesh,
394  nx_divs, ny_divs, nz_divs,
395  length,
396  height,
397  width,
398  e_type);
399 
400  // we now uniformly refine this mesh
401  if (n_refine) {
402 
403  libMesh::MeshRefinement(mesh).uniformly_refine(n_refine);
404  libMesh::MeshTools::Modification::flatten(mesh);
405  }
406  }
407 
408 
409 
410  template <typename Context>
411  inline void
412  init_analysis_dirichlet_conditions(Context& c) {
413 
414  c.sys->get_dof_map().add_dirichlet_boundary
415  (libMesh::DirichletBoundary({1}, {0, 1, 2}, libMesh::ZeroFunction<real_t>()));
416  }
417 
418 
419 
420  template <typename ScalarType, typename InitType>
421  std::unique_ptr<pressure_t<ScalarType>>
422  build_pressure_load(InitType& c) {
423 
424  real_t
425  length = c.input("length", "length of domain along x-axis", 0.3),
426  frac = c.input("loadlength_fraction", "fraction of boundary length on which pressure will act", 0.125),
427  p_val = c.input("pressure", "pressure on side of domain", 5.e7);
428  c.p_side_id = 7;
429 
430  std::unique_ptr<pressure_t<ScalarType>>
431  press(new pressure_t<ScalarType>(p_val, length, frac));
432 
433  return press;
434  }
435 
436 
437 
438  template <typename ScalarType, typename Context>
439  inline void
440  init_simp_dvs
441  (Context &c,
442  MAST::Optimization::DesignParameterVector<ScalarType> &dvs) {
443 
444  //
445  // this assumes that density variable has a constant value per element
446  //
447  Assert2(c.rho_fe_family == libMesh::LAGRANGE,
448  c.rho_fe_family, libMesh::LAGRANGE,
449  "Method assumes Lagrange interpolation function for density");
450  Assert2(c.rho_fe_order == libMesh::FIRST,
451  c.rho_fe_order, libMesh::FIRST,
452  "Method assumes Lagrange interpolation function for density");
453 
454  real_t
455  tol = 1.e-12,
456  rho_min = c.input("rho_min", "lower limit on density variable", 0.),
457  vf = c.input("volume_fraction", "upper limit for the volume fraction", 0.2);
458 
459  uint_t
460  sys_num = c.rho_sys->number(),
461  first_dof = c.rho_sys->get_dof_map().first_dof(c.rho_sys->comm().rank()),
462  end_dof = c.rho_sys->get_dof_map().end_dof(c.rho_sys->comm().rank()),
463  dof_id = 0;
464 
465  real_t
466  val = 0.;
467 
468  libMesh::MeshBase::const_element_iterator
469  e_it = c.mesh->active_local_elements_begin(),
470  e_end = c.mesh->active_local_elements_end();
471 
472  std::set<const libMesh::Node*> nodes;
473 
474  for ( ; e_it != e_end; e_it++) {
475 
476  const libMesh::Elem* e = *e_it;
477 
478  Assert0(e->type() == libMesh::HEX8 ||
479  e->type() == libMesh::HEX27,
480  "Method requires Hex8/Hex27 element");
481 
482  // only the first eight nodes of the hex element are used
483  for (uint_t i=0; i<8; i++) {
484 
485  const libMesh::Node& n = *e->node_ptr(i);
486 
487  // if we have alredy operated on this node, then
488  // we skip it
489  if (nodes.count(&n))
490  continue;
491 
492  // otherwise, we add it to the set of operated nodes and
493  // check if a design parameter should be computed for this
494  nodes.insert(&n);
495 
496  dof_id = n.dof_number(sys_num, 0, 0);
497 
498  MAST::Optimization::DesignParameter<ScalarType>
499  *dv = new MAST::Optimization::DesignParameter<ScalarType>(vf);
500  dv->set_point(n(0), n(1), n(2));
501 
502  if (dof_id >= first_dof &&
503  dof_id < end_dof)
504  dvs.add_topology_parameter(*dv, dof_id);
505  else
506  dvs.add_ghosted_topology_parameter(*dv, dof_id);
507  }
508  }
509 
510  dvs.synchronize(c.rho_sys->get_dof_map());
511  c.rho_sys->solution->close();
512  }
513 };
514 
515 } // namespace Generation
516 } // namespace Mesh
517 } // namespace MAST
518 
519 
520 #endif // __mast_mesh_generation_bracket_3d_h__
Error
#define Error(cond, msg)
Definition: exceptions.hpp:166
MAST::Mesh::Generation::Bracket3D::init_simp_dvs
void init_simp_dvs(Context &c, MAST::Optimization::DesignParameterVector< ScalarType > &dvs)
Definition: bracket3d.hpp:441
MAST::Optimization::DesignParameterVector::add_ghosted_topology_parameter
MAST::Base::ParameterData & add_ghosted_topology_parameter(MAST::Optimization::DesignParameter< ScalarType > &p, const uint_t id)
Definition: design_parameter_vector.hpp:212
MAST::Mesh::Generation::Bracket3D::idx
uint_t idx(const libMesh::ElemType type, const uint_t nx, const uint_t ny, const uint_t i, const uint_t j, const uint_t k)
Definition: bracket3d.hpp:94
MAST::Mesh::Generation::Bracket3D::init_analysis_dirichlet_conditions
void init_analysis_dirichlet_conditions(Context &c)
Definition: bracket3d.hpp:412
MAST::Mesh::Generation::Bracket3D::Pressure::value
ScalarType value(ContextType &c) const
Definition: bracket3d.hpp:62
MAST::Mesh::Generation::Bracket3D::init_analysis_mesh
void init_analysis_mesh(Context &c, libMesh::UnstructuredMesh &mesh)
Definition: bracket3d.hpp:358
MAST::Optimization::DesignParameterVector::add_topology_parameter
MAST::Base::ParameterData & add_topology_parameter(MAST::Optimization::DesignParameter< ScalarType > &p, const uint_t id)
Definition: design_parameter_vector.hpp:189
parameter_data.hpp
MAST::Mesh::Generation::Bracket3D::Pressure::derivative
ScalarType derivative(ContextType &c, const ScalarFieldType &f) const
Definition: bracket3d.hpp:68
mast_data_types.h
getpot_wrapper.hpp
Assert0
#define Assert0(cond, msg)
Definition: exceptions.hpp:134
MAST::Mesh::Generation::Bracket3D::Pressure::Pressure
Pressure(real_t p, real_t l1, real_t frac)
Definition: bracket3d.hpp:54
uint_t
unsigned int uint_t
Definition: mast_data_types.h:43
Assert2
#define Assert2(cond, v1, v2, msg)
Definition: exceptions.hpp:152
design_parameter_vector.hpp
MAST::Optimization::DesignParameter::set_point
void set_point(real_t x, real_t y=0., real_t z=0.)
Definition: design_parameter.hpp:49
MAST::Mesh::Generation::Bracket3D::build_pressure_load
std::unique_ptr< pressure_t< ScalarType > > build_pressure_load(InitType &c)
Definition: bracket3d.hpp:422
real_t
double real_t
Definition: mast_data_types.h:42
MAST::Mesh::Generation::Bracket3D::build_cube
void build_cube(libMesh::UnstructuredMesh &mesh, const uint_t nx, const uint_t ny, const uint_t nz, const real_t length, const real_t height, const real_t width, const libMesh::ElemType type)
Definition: bracket3d.hpp:121
MAST::Optimization::DesignParameterVector::synchronize
void synchronize(const libMesh::DofMap &dof_map)
Definition: design_parameter_vector.hpp:262
MAST::Numerics::Utility::build
std::unique_ptr< ValType > build(const libMesh::System &sys)
Definition: utility.hpp:284
MAST::Mesh::Generation::Bracket3D::reference_volume
real_t reference_volume(Context &c)
Definition: bracket3d.hpp:83