MAST3
Multidisciplinary-design Adaptation and Sensitivity Toolkit (MAST)
panel3d.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_panel_3d_h__
21 #define __mast_mesh_generation_panel_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 Panel3D {
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=(fabs(c.qp_location(0)-_l1*0.5) <= 0.5*_frac*_l1)?_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::Truss3D::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.03),
88  width = c.input("width", "length of domain along z-axis", 0.3);
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 real_t dirichlet_length_fraction,
129  const libMesh::ElemType type) {
130 
131  Assert0(type == libMesh::HEX8 || type == libMesh::HEX27,
132  "Method only implemented for Hex8/Hex27");
133 
134  // Clear the mesh and start from scratch
135  mesh.clear();
136 
137  libMesh::BoundaryInfo & boundary_info = mesh.get_boundary_info();
138 
139  mesh.set_mesh_dimension(3);
140  mesh.set_spatial_dimension(3);
141  mesh.reserve_elem(nx*ny*nz);
142 
143  if (type == libMesh::HEX8)
144  mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
145  else if (type == libMesh::HEX27)
146  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
147 
148  real_t
149  xmax = length,
150  ymax = height,
151  zmax = width;
152 
153 
154 
155  std::map<uint_t, libMesh::Node*> nodes;
156 
157  // Build the nodes.
158  uint_t
159  node_id = 0,
160  n = 0;
161  switch (type)
162  {
163  case libMesh::HEX8: {
164 
165  for (uint_t k=0; k<=nz; k++)
166  for (uint_t j=0; j<=ny; j++)
167  for (uint_t i=0; i<=nx; i++) {
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  node_id++;
174  }
175 
176 
177  break;
178  }
179 
180  case libMesh::HEX27: {
181 
182  for (uint_t k=0; k<=(2*nz); k++)
183  for (uint_t j=0; j<=(2*ny); j++)
184  for (uint_t i=0; i<=(2*nx); i++) {
185  nodes[node_id] =
186  mesh.add_point(libMesh::Point(static_cast<real_t>(i)/static_cast<real_t>(2*nx)*length,
187  static_cast<real_t>(j)/static_cast<real_t>(2*ny)*height,
188  static_cast<real_t>(k)/static_cast<real_t>(2*nz)*width),
189  n++);
190  node_id++;
191  }
192 
193  break;
194  }
195 
196 
197  default:
198  Assert0(false, "ERROR: Unrecognized 3D element type.");
199  }
200 
201  // Build the elements.
202  uint_t
203  elem_id = 0;
204  switch (type)
205  {
206  case libMesh::HEX8:
207  {
208  for (uint_t k=0; k<nz; k++)
209  for (uint_t j=0; j<ny; j++)
210  for (uint_t i=0; i<nx; i++) {
211 
212  libMesh::Elem
213  *elem = libMesh::Elem::build(libMesh::HEX8).release();
214  elem->set_id(elem_id++);
215  mesh.add_elem(elem);
216 
217  elem->set_node(0) = nodes[idx(type,nx,ny,i,j,k) ];
218  elem->set_node(1) = nodes[idx(type,nx,ny,i+1,j,k) ];
219  elem->set_node(2) = nodes[idx(type,nx,ny,i+1,j+1,k) ];
220  elem->set_node(3) = nodes[idx(type,nx,ny,i,j+1,k) ];
221  elem->set_node(4) = nodes[idx(type,nx,ny,i,j,k+1) ];
222  elem->set_node(5) = nodes[idx(type,nx,ny,i+1,j,k+1) ];
223  elem->set_node(6) = nodes[idx(type,nx,ny,i+1,j+1,k+1)];
224  elem->set_node(7) = nodes[idx(type,nx,ny,i,j+1,k+1) ];
225 
226  if (k == 0)
227  boundary_info.add_side(elem, 0, 0);
228 
229  if (k == (nz-1))
230  boundary_info.add_side(elem, 5, 5);
231 
232  if (j == 0)
233  boundary_info.add_side(elem, 1, 1);
234 
235  if (j == (ny-1))
236  boundary_info.add_side(elem, 3, 3);
237 
238  if (i == 0)
239  boundary_info.add_side(elem, 4, 4);
240 
241  if (i == (nx-1))
242  boundary_info.add_side(elem, 2, 2);
243 
244  // dirichlet on face 0
245  if (k == 0 && j >= (1.-dirichlet_length_fraction)* ny)
246  boundary_info.add_side(elem, 0, 6);
247 
248  // dirichlet on face 2
249  if (i == (nx-1) && j >= (1.-dirichlet_length_fraction)* ny)
250  boundary_info.add_side(elem, 2, 7);
251 
252  // dirichlet on face 4
253  if (i == 0 && j >= (1.-dirichlet_length_fraction)* ny)
254  boundary_info.add_side(elem, 4, 8);
255 
256  // dirichlet on face 5
257  if (k == (nz-1) && j >= (1.-dirichlet_length_fraction)* ny)
258  boundary_info.add_side(elem, 5, 9);
259  }
260  break;
261  }
262 
263 
264  case libMesh::HEX27: {
265 
266  for (uint_t k=0; k<(2*nz); k += 2)
267  for (uint_t j=0; j<(2*ny); j += 2)
268  for (uint_t i=0; i<(2*nx); i += 2)
269  {
270  libMesh::Elem
271  *elem = libMesh::Elem::build(libMesh::HEX27).release();
272  elem->set_id(elem_id++);
273  mesh.add_elem(elem);
274 
275  elem->set_node(0) = nodes[idx(type,nx,ny,i, j, k) ];
276  elem->set_node(1) = nodes[idx(type,nx,ny,i+2,j, k) ];
277  elem->set_node(2) = nodes[idx(type,nx,ny,i+2,j+2,k) ];
278  elem->set_node(3) = nodes[idx(type,nx,ny,i, j+2,k) ];
279  elem->set_node(4) = nodes[idx(type,nx,ny,i, j, k+2)];
280  elem->set_node(5) = nodes[idx(type,nx,ny,i+2,j, k+2)];
281  elem->set_node(6) = nodes[idx(type,nx,ny,i+2,j+2,k+2)];
282  elem->set_node(7) = nodes[idx(type,nx,ny,i, j+2,k+2)];
283  elem->set_node(8) = nodes[idx(type,nx,ny,i+1,j, k) ];
284  elem->set_node(9) = nodes[idx(type,nx,ny,i+2,j+1,k) ];
285  elem->set_node(10) = nodes[idx(type,nx,ny,i+1,j+2,k) ];
286  elem->set_node(11) = nodes[idx(type,nx,ny,i, j+1,k) ];
287  elem->set_node(12) = nodes[idx(type,nx,ny,i, j, k+1)];
288  elem->set_node(13) = nodes[idx(type,nx,ny,i+2,j, k+1)];
289  elem->set_node(14) = nodes[idx(type,nx,ny,i+2,j+2,k+1)];
290  elem->set_node(15) = nodes[idx(type,nx,ny,i, j+2,k+1)];
291  elem->set_node(16) = nodes[idx(type,nx,ny,i+1,j, k+2)];
292  elem->set_node(17) = nodes[idx(type,nx,ny,i+2,j+1,k+2)];
293  elem->set_node(18) = nodes[idx(type,nx,ny,i+1,j+2,k+2)];
294  elem->set_node(19) = nodes[idx(type,nx,ny,i, j+1,k+2)];
295 
296  elem->set_node(20) = nodes[idx(type,nx,ny,i+1,j+1,k) ];
297  elem->set_node(21) = nodes[idx(type,nx,ny,i+1,j, k+1)];
298  elem->set_node(22) = nodes[idx(type,nx,ny,i+2,j+1,k+1)];
299  elem->set_node(23) = nodes[idx(type,nx,ny,i+1,j+2,k+1)];
300  elem->set_node(24) = nodes[idx(type,nx,ny,i, j+1,k+1)];
301  elem->set_node(25) = nodes[idx(type,nx,ny,i+1,j+1,k+2)];
302  elem->set_node(26) = nodes[idx(type,nx,ny,i+1,j+1,k+1)];
303 
304  if (k == 0)
305  boundary_info.add_side(elem, 0, 0);
306 
307  if (k == 2*(nz-1))
308  boundary_info.add_side(elem, 5, 5);
309 
310  if (j == 0)
311  boundary_info.add_side(elem, 1, 1);
312 
313  if (j == 2*(ny-1))
314  boundary_info.add_side(elem, 3, 3);
315 
316  if (i == 0)
317  boundary_info.add_side(elem, 4, 4);
318 
319  if (i == 2*(nx-1))
320  boundary_info.add_side(elem, 2, 2);
321 
322  // dirichlet on face 0
323  if (k == 0 && j >= (1.-dirichlet_length_fraction)* 2*ny)
324  boundary_info.add_side(elem, 0, 6);
325 
326  // dirichlet on face 2
327  if (i == 2*(nx-1) && j >= (1.-dirichlet_length_fraction)* 2*ny)
328  boundary_info.add_side(elem, 2, 7);
329 
330  // dirichlet on face 4
331  if (i == 0 && j >= (1.-dirichlet_length_fraction)* 2*ny)
332  boundary_info.add_side(elem, 4, 8);
333 
334  // dirichlet on face 5
335  if (k == 2*(nz-1) && j >= (1.-dirichlet_length_fraction)* 2*ny)
336  boundary_info.add_side(elem, 5, 9);
337 
338  }
339  break;
340  }
341 
342  default:
343  Assert0(false, "ERROR: Unrecognized 3D element type.");
344  }
345 
346  // Add sideset names to boundary info (Z axis out of the screen)
347  boundary_info.sideset_name(0) = "back";
348  boundary_info.sideset_name(1) = "bottom";
349  boundary_info.sideset_name(2) = "right";
350  boundary_info.sideset_name(3) = "top";
351  boundary_info.sideset_name(4) = "left";
352  boundary_info.sideset_name(5) = "front";
353  boundary_info.sideset_name(6) = "back_dirichlet";
354  boundary_info.sideset_name(7) = "right_dirichlet";
355  boundary_info.sideset_name(8) = "left_dirichlet";
356  boundary_info.sideset_name(9) = "front_dirichlet";
357 
358  // Add nodeset names to boundary info
359  boundary_info.nodeset_name(0) = "back";
360  boundary_info.nodeset_name(1) = "bottom";
361  boundary_info.nodeset_name(2) = "right";
362  boundary_info.nodeset_name(3) = "top";
363  boundary_info.nodeset_name(4) = "left";
364  boundary_info.nodeset_name(5) = "front";
365 
366  // Done building the mesh. Now prepare it for use.
367  mesh.prepare_for_use ();
368  }
369 
370 
371 
372  template <typename Context>
373  inline void
374  init_analysis_mesh(Context& c,
375  libMesh::UnstructuredMesh& mesh) {
376 
377  real_t
378  length = c.input("length", "length of domain along x-axis", 0.3),
379  height = c.input("height", "length of domain along y-axis", 0.03),
380  width = c.input("width", "length of domain along z-axis", 0.3),
381  dirichlet_length_fraction = c.input
382  ("dirichlet_length_fraction",
383  "length fraction of the truss boundary where dirichlet condition is applied",
384  0.1);
385 
386  uint_t
387  nx_divs = c.input("nx_divs", "number of elements along x-axis", 100),
388  ny_divs = c.input("ny_divs", "number of elements along y-axis", 10),
389  nz_divs = c.input("nz_divs", "number of elements along z-axis", 100),
390  n_refine= c.input("n_uniform_refinement", "number of times the mesh is uniformly refined", 0);
391 
392  std::string
393  t = c.input("elem_type", "type of geometric element in the mesh", "hex8");
394 
395  libMesh::ElemType
396  e_type = libMesh::Utility::string_to_enum<libMesh::ElemType>(t);
397 
398  //
399  // if high order FE is used, libMesh requires atleast a second order
400  // geometric element.
401  //
402  if (c.sol_fe_order > 1 && e_type == libMesh::HEX8)
403  e_type = libMesh::HEX27;
404  else if (c.rho_fe_order > 1 && e_type == libMesh::TET4)
405  e_type = libMesh::TET10;
406 
407  //
408  // initialize the mesh with one element
409  //
410  build_cube(mesh,
411  nx_divs, ny_divs, nz_divs,
412  length,
413  height,
414  width,
415  dirichlet_length_fraction,
416  e_type);
417 
418  // we now uniformly refine this mesh
419  if (n_refine) {
420 
421  libMesh::MeshRefinement(mesh).uniformly_refine(n_refine);
422  libMesh::MeshTools::Modification::flatten(mesh);
423  }
424  }
425 
426 
427 
428  template <typename Context>
429  inline void
430  init_analysis_dirichlet_conditions(Context& c) {
431 
432  c.sys->get_dof_map().add_dirichlet_boundary
433  (libMesh::DirichletBoundary({6, 7, 8, 9}, {0, 1, 2}, libMesh::ZeroFunction<real_t>()));
434  }
435 
436 
437 
438  template <typename ScalarType, typename InitType>
439  std::unique_ptr<pressure_t<ScalarType>>
440  build_pressure_load(InitType& c) {
441 
442  real_t
443  length = c.input("length", "length of domain along x-axis", 0.3),
444  frac = c.input("loadlength_fraction", "fraction of boundary length on which pressure will act", 1.0),
445  p_val = c.input("pressure", "pressure on side of domain", 2.e6);
446  c.p_side_id = 3;
447 
448  std::unique_ptr<pressure_t<ScalarType>>
449  press(new pressure_t<ScalarType>(p_val, length, frac));
450 
451  return press;
452  }
453 
454 
455 
456  template <typename ScalarType, typename Context>
457  inline void
458  init_simp_dvs
459  (Context &c,
460  MAST::Optimization::DesignParameterVector<ScalarType> &dvs) {
461 
462  //
463  // this assumes that density variable has a constant value per element
464  //
465  Assert2(c.rho_fe_family == libMesh::LAGRANGE,
466  c.rho_fe_family, libMesh::LAGRANGE,
467  "Method assumes Lagrange interpolation function for density");
468  Assert2(c.rho_fe_order == libMesh::FIRST,
469  c.rho_fe_order, libMesh::FIRST,
470  "Method assumes Lagrange interpolation function for density");
471 
472  real_t
473  vf = c.input("volume_fraction", "upper limit for the volume fraction", 0.2);
474 
475  uint_t
476  sys_num = c.rho_sys->number(),
477  first_dof = c.rho_sys->get_dof_map().first_dof(c.rho_sys->comm().rank()),
478  end_dof = c.rho_sys->get_dof_map().end_dof(c.rho_sys->comm().rank()),
479  dof_id = 0;
480 
481  real_t
482  val = 0.;
483 
484  libMesh::MeshBase::const_element_iterator
485  e_it = c.mesh->active_local_elements_begin(),
486  e_end = c.mesh->active_local_elements_end();
487 
488  std::set<const libMesh::Node*> nodes;
489 
490  for ( ; e_it != e_end; e_it++) {
491 
492  const libMesh::Elem* e = *e_it;
493 
494  Assert0(e->type() == libMesh::HEX8 ||
495  e->type() == libMesh::HEX27,
496  "Method requires Hex8/Hex27 element");
497 
498  // only the first eight nodes of the hex element are used
499  for (uint_t i=0; i<8; i++) {
500 
501  const libMesh::Node& n = *e->node_ptr(i);
502 
503  // if we have alredy operated on this node, then
504  // we skip it
505  if (nodes.count(&n))
506  continue;
507 
508  // otherwise, we add it to the set of operated nodes and
509  // check if a design parameter should be computed for this
510  nodes.insert(&n);
511 
512  dof_id = n.dof_number(sys_num, 0, 0);
513 
514  MAST::Optimization::DesignParameter<ScalarType>
515  *dv = new MAST::Optimization::DesignParameter<ScalarType>(vf);
516  dv->set_point(n(0), n(1), n(2));
517 
518  if (dof_id >= first_dof &&
519  dof_id < end_dof)
520  dvs.add_topology_parameter(*dv, dof_id);
521  else
522  dvs.add_ghosted_topology_parameter(*dv, dof_id);
523  }
524  }
525 
526  dvs.synchronize(c.rho_sys->get_dof_map());
527  c.rho_sys->solution->close();
528  }
529 };
530 
531 } // namespace Generation
532 } // namespace Mesh
533 } // namespace MAST
534 
535 
536 #endif // __mast_mesh_generation_truss_3d_h__
Error
#define Error(cond, msg)
Definition: exceptions.hpp:166
MAST::Mesh::Generation::Panel3D::Pressure::value
ScalarType value(ContextType &c) const
Definition: panel3d.hpp:62
MAST::Mesh::Generation::Panel3D::init_analysis_dirichlet_conditions
void init_analysis_dirichlet_conditions(Context &c)
Definition: panel3d.hpp:430
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::Panel3D::build_pressure_load
std::unique_ptr< pressure_t< ScalarType > > build_pressure_load(InitType &c)
Definition: panel3d.hpp:440
MAST::Mesh::Generation::Panel3D::reference_volume
real_t reference_volume(Context &c)
Definition: panel3d.hpp:83
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
MAST::Mesh::Generation::Panel3D::init_analysis_mesh
void init_analysis_mesh(Context &c, libMesh::UnstructuredMesh &mesh)
Definition: panel3d.hpp:374
parameter_data.hpp
MAST::Mesh::Generation::Panel3D::dim
static const uint_t dim
Definition: panel3d.hpp:77
mast_data_types.h
MAST::Mesh::Generation::Panel3D::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: panel3d.hpp:94
getpot_wrapper.hpp
MAST::Mesh::Generation::Panel3D::init_simp_dvs
void init_simp_dvs(Context &c, MAST::Optimization::DesignParameterVector< ScalarType > &dvs)
Definition: panel3d.hpp:459
Assert0
#define Assert0(cond, msg)
Definition: exceptions.hpp:134
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
real_t
double real_t
Definition: mast_data_types.h:42
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::Panel3D::Pressure::derivative
ScalarType derivative(ContextType &c, const ScalarFieldType &f) const
Definition: panel3d.hpp:68
MAST::Mesh::Generation::Panel3D::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 real_t dirichlet_length_fraction, const libMesh::ElemType type)
Definition: panel3d.hpp:121
MAST::Mesh::Generation::Panel3D::Pressure::Pressure
Pressure(real_t p, real_t l1, real_t frac)
Definition: panel3d.hpp:54