diff --git a/examples/structural/base/bracket_2d_model.h b/examples/structural/base/bracket_2d_model.h index 57f39382..35e0af12 100644 --- a/examples/structural/base/bracket_2d_model.h +++ b/examples/structural/base/bracket_2d_model.h @@ -49,10 +49,13 @@ class DisciplineBase; class BoundaryConditionBase; class FunctionBase; class Parameter; +class Solid2DSectionElementPropertyCard; namespace Examples { struct Bracket2DModel { + using SectionPropertyCardType = MAST::Solid2DSectionElementPropertyCard; + template static Real reference_volume(Opt& opt); @@ -637,7 +640,7 @@ MAST::Examples::Bracket2DModel::initialize_level_set_solution(Opt& opt) { nx_m = opt._input("level_set_nx_divs", "number of elements of level-set mesh along x-axis", 10), ny_m = opt._input("level_set_ny_divs", "number of elements of level-set mesh along y-axis", 10); - MAST::Examples::LevelSetNucleationFunction + MAST::Examples::LevelSetNucleationFunction2D phi(0., 0., length, height, nx_m, ny_m, nx_h, ny_h); opt._level_set_sys_init->initialize_solution(phi); diff --git a/examples/structural/base/bracket_3d_model.h b/examples/structural/base/bracket_3d_model.h new file mode 100644 index 00000000..f2ea3ff1 --- /dev/null +++ b/examples/structural/base/bracket_3d_model.h @@ -0,0 +1,660 @@ +/* + * MAST: Multidisciplinary-design Adaptation and Sensitivity Toolkit + * Copyright (C) 2013-2020 Manav Bhatia and MAST authors + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * This library is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with this library; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + */ + +#ifndef __mast_topology_3d_bracket_model__ +#define __mast_topology_3d_bracket_model__ + +// MAST includes +#include "base/mast_data_types.h" +#include "examples/base/input_wrapper.h" +#include "examples/structural/base/level_set_nucleation.h" +#include "base/boundary_condition_base.h" +#include "base/field_function_base.h" +#include "base/physics_discipline_base.h" +#include "boundary_condition/dirichlet_boundary_condition.h" +#include "level_set/level_set_parameter.h" + +// libMesh includes +#include "libmesh/system.h" +#include "libmesh/unstructured_mesh.h" +#include "libmesh/fe_type.h" +#include "libmesh/string_to_enum.h" +#include "libmesh/mesh_generation.h" +#include "libmesh/elem.h" +#include "libmesh/node.h" + + + +namespace MAST { + +// Forward declerations +class DisciplineBase; +class BoundaryConditionBase; +class FunctionBase; +class Parameter; +class IsotropicElementPropertyCard3D; + +namespace Examples { +struct Bracket3DModel { + + using SectionPropertyCardType = MAST::IsotropicElementPropertyCard3D; + + template + static Real reference_volume(Opt& opt); + + template + static void init_analysis_mesh(Opt& opt, libMesh::UnstructuredMesh& mesh); + + template + static void init_level_set_mesh(Opt& opt, libMesh::UnstructuredMesh& mesh); + + template + static void init_analysis_dirichlet_conditions(Opt& opt); + + template + static void init_indicator_dirichlet_conditions(Opt& opt); + + template + static void init_structural_loads(Opt& opt); + + template + static MAST::BoundaryConditionBase& + init_structural_shifted_boudnary_load(Opt& opt, unsigned int bid); + + template + static void init_indicator_loads(Opt& opt); + + template + static void init_level_set_dvs(Opt& opt); + + template + static void initialize_level_set_solution(Opt& opt); + + template + static void init_simp_dvs(Opt& opt); + + template + static void _delete_elems_from_bracket_mesh(Opt& opt, libMesh::MeshBase &mesh); + + class BracketLoad: + public MAST::FieldFunction { + public: + BracketLoad(const std::string& nm, Real p, Real l1, Real fraction): + MAST::FieldFunction(nm), _p(p), _l1(l1), _frac(fraction) { } + ~BracketLoad() {} + void operator() (const libMesh::Point& p, const Real t, Real& v) const { + if (fabs(p(0) >= _l1*(1.-_frac))) v = _p; + else v = 0.; + } + void derivative(const MAST::FunctionBase& f, const libMesh::Point& p, const Real t, Real& v) const { + v = 0.; + } + protected: + Real _p, _l1, _frac; + }; +}; + +} +} + + +template +Real +MAST::Examples::Bracket3DModel:: +reference_volume(Opt& opt) { + + Real + length = opt._input("length", "length of domain along x-axis", 0.3), + height = opt._input("height", "length of domain along y-axis", 0.3), + width = opt._input("width", "length of domain along z-axis", 0.3); + + return length * height * width; +} + + + +template +void +MAST::Examples::Bracket3DModel:: +init_analysis_mesh(Opt& opt, + libMesh::UnstructuredMesh& mesh) { + + Real + length = opt._input("length", "length of domain along x-axis", 0.3), + height = opt._input("height", "length of domain along y-axis", 0.3), + width = opt._input("width", "length of domain along z-axis", 0.3); + + unsigned int + nx_divs = opt._input("nx_divs", "number of elements along x-axis", 20), + ny_divs = opt._input("ny_divs", "number of elements along y-axis", 20), + nz_divs = opt._input("nz_divs", "number of elements along z-axis", 20); + + if (nx_divs%10 != 0 || ny_divs%10 != 0) libmesh_error(); + + std::string + t = opt._input("elem_type", "type of geometric element in the mesh", "hex8"); + + libMesh::ElemType + e_type = libMesh::Utility::string_to_enum(t); + + // + // if high order FE is used, libMesh requires atleast a second order + // geometric element. + // + if (opt._fetype.order > 1 && e_type == libMesh::HEX8) + e_type = libMesh::HEX27; + else if (opt._fetype.order > 1 && e_type == libMesh::TET4) + e_type = libMesh::TET10; + + // + // initialize the mesh with one element + // + libMesh::MeshTools::Generation::build_cube(mesh, + nx_divs, ny_divs, nz_divs, + 0, length, + 0, height, + 0, width, + e_type); + + _delete_elems_from_bracket_mesh(opt, mesh); +} + + +template +void +MAST::Examples::Bracket3DModel:: +init_level_set_mesh(Opt& opt, + libMesh::UnstructuredMesh& mesh) { + + Real + length = opt._input("length", "length of domain along x-axis", 0.3), + height = opt._input("height", "length of domain along y-axis", 0.3), + width = opt._input("width", "length of domain along z-axis", 0.3); + + unsigned int + nx_divs = opt._input("level_set_nx_divs", "number of elements of level-set mesh along x-axis", 10), + ny_divs = opt._input("level_set_ny_divs", "number of elements of level-set mesh along y-axis", 10), + nz_divs = opt._input("level_set_nz_divs", "number of elements of level-set mesh along z-axis", 10); + + if (nx_divs%10 != 0 || ny_divs%10 != 0) libmesh_error(); + + libMesh::ElemType + e_type = libMesh::HEX8; + + // initialize the mesh with one element + libMesh::MeshTools::Generation::build_cube(mesh, + nx_divs, ny_divs, nz_divs, + 0, length, + 0, height, + 0, width, + e_type); + + _delete_elems_from_bracket_mesh(opt, mesh); +} + + + +template +void +MAST::Examples::Bracket3DModel:: +init_analysis_dirichlet_conditions(Opt& opt) { + + MAST::DirichletBoundaryCondition + *dirichlet = new MAST::DirichletBoundaryCondition; // bottom boundary + dirichlet->init(1, opt._sys_init->vars()); + opt._discipline->add_dirichlet_bc(1, *dirichlet); + opt._boundary_conditions.insert(dirichlet); + + opt._discipline->init_system_dirichlet_bc(*opt._sys); +} + + + +template +void +MAST::Examples::Bracket3DModel:: +init_indicator_dirichlet_conditions(Opt& opt) { + + MAST::DirichletBoundaryCondition + *dirichlet = new MAST::DirichletBoundaryCondition; // bottom boundary + dirichlet->init(1, opt._indicator_sys_init->vars()); + opt._indicator_discipline->add_dirichlet_bc(1, *dirichlet); + opt._boundary_conditions.insert(dirichlet); + + opt._indicator_discipline->init_system_dirichlet_bc(*opt._indicator_sys); + opt._dirichlet_bc_ids.insert(1); +} + + + +template +void +MAST::Examples::Bracket3DModel::init_structural_loads(Opt& opt) { + + Real + length = opt._input("length", "length of domain along x-axis", 0.3), + frac = opt._input("loadlength_fraction", "fraction of boundary length on which pressure will act", 0.125), + p_val = opt._input("pressure", "pressure on side of domain", 5.e7); + + BracketLoad + *press_f = new BracketLoad( "pressure", p_val, length, frac); + + // + // initialize the load + // + MAST::BoundaryConditionBase + *p_load = new MAST::BoundaryConditionBase(MAST::SURFACE_PRESSURE); + + p_load->add(*press_f); + opt._discipline->add_side_load(7, *p_load); + opt._boundary_conditions.insert(p_load); + + opt._field_functions.insert(press_f); +} + + +template +MAST::BoundaryConditionBase& +MAST::Examples::Bracket3DModel::init_structural_shifted_boudnary_load(Opt& opt, + unsigned int bid) { + + class ZeroTraction: public MAST::FieldFunction { + public: + ZeroTraction(): MAST::FieldFunction("traction") {} + virtual ~ZeroTraction() {} + virtual void operator() (const libMesh::Point& pt, const Real t, RealVectorX& v) const {v.setZero(3);} + virtual void derivative(const MAST::FunctionBase& f, const libMesh::Point& pt, const Real t, RealVectorX& v) const + {v.setZero(3);} + }; + + ZeroTraction + *trac_f = new ZeroTraction; + + MAST::BoundaryConditionBase + *load = new MAST::BoundaryConditionBase(MAST::SURFACE_TRACTION_SHIFTED_BOUNDARY); + + load->add(*opt._level_set_vel); + load->add(*trac_f); + opt._discipline->add_side_load(bid, *load); + opt._boundary_conditions.insert(load); + + opt._field_functions.insert(trac_f); + return *load; +} + + +template +void +MAST::Examples::Bracket3DModel::init_indicator_loads(Opt& opt) { + + Real + length = opt._input("length", "length of domain along x-axis", 0.3), + frac = opt._input("loadlength_fraction", "fraction of boundary length on which pressure will act", 0.125); + + BracketLoad + *flux_f = new BracketLoad("heat_flux", -2.e6, length, frac); + + // + // initialize the load + // + MAST::BoundaryConditionBase + *f_load = new MAST::BoundaryConditionBase(MAST::HEAT_FLUX); + + f_load->add(*flux_f); + opt._indicator_discipline->add_side_load(7, *f_load); + opt._boundary_conditions.insert(f_load); + + opt._field_functions.insert(flux_f); +} + + + +template +void +MAST::Examples::Bracket3DModel::init_level_set_dvs(Opt& opt) { + + libmesh_assert(opt._initialized); + // + // this assumes that level set is defined using lagrange shape functions + // + libmesh_assert_equal_to(opt._level_set_fetype.family, libMesh::LAGRANGE); + + Real + tol = 1.e-12, + l_frac = 0.4,//_input("length_fraction", "fraction of length along x-axis that is in the bracket", 0.4), + h_frac = 0.4,//_input( "height_fraction", "fraction of length along y-axis that is in the bracket", 0.4), + length = opt._input("length", "length of domain along x-axis", 0.3), + height = opt._input("height", "length of domain along y-axis", 0.3), + x_lim = length * l_frac, + y_lim = height * (1.-h_frac), + frac = opt._input("loadlength_fraction", "fraction of boundary length on which pressure will act", 0.125), + filter_radius = opt._input("filter_radius", "radius of geometric filter for level set field", 0.015); + + unsigned int + dof_id = 0, + n_vars = 0; + + Real + val = 0.; + + // + // all ranks will have DVs defined for all variables. So, we should be + // operating on a replicated mesh + // + libmesh_assert(opt._level_set_mesh->is_replicated()); + + std::vector local_phi(opt._level_set_sys->solution->size()); + opt._level_set_sys->solution->localize(local_phi); + + // + // iterate over all the node values + // + libMesh::MeshBase::const_node_iterator + it = opt._level_set_mesh->nodes_begin(), + end = opt._level_set_mesh->nodes_end(); + + // + // maximum number of dvs is the number of nodes on the level set function + // mesh. We will evaluate the actual number of dvs + // + opt._dv_params.reserve(opt._level_set_mesh->n_nodes()); + n_vars = 0; + + for ( ; it!=end; it++) { + + const libMesh::Node& n = **it; + + dof_id = n.dof_number(0, 0, 0); + + if ((n(1)-filter_radius) <= y_lim && + (n(0)+filter_radius) >= length*(1.-frac)) { + + // + // set value at the constrained points to a small positive number + // material here + // + if (dof_id >= opt._level_set_sys->solution->first_local_index() && + dof_id < opt._level_set_sys->solution->last_local_index()) + opt._level_set_sys->solution->set(dof_id, 1.e0); + } + else { + + std::ostringstream oss; + oss << "dv_" << n_vars; + val = local_phi[dof_id]; + + // + // on the boundary, set everything to be zero, so that there + // is always a boundary there that the optimizer can move + // + if (n(0) < tol || // left boundary + std::fabs(n(0) - length) < tol || // right boundary + std::fabs(n(1) - height) < tol || // top boundary + (n(0) >= x_lim && n(1) <= y_lim)) { + + if (dof_id >= opt._level_set_sys->solution->first_local_index() && + dof_id < opt._level_set_sys->solution->last_local_index()) + opt._level_set_sys->solution->set(dof_id, -1.0); + val = -1.0; + } + + opt._dv_params.push_back(std::pair()); + opt._dv_params[n_vars].first = dof_id; + opt._dv_params[n_vars].second = new MAST::LevelSetParameter(oss.str(), val, &n); + opt._dv_params[n_vars].second->set_as_topology_parameter(true); + opt._dv_dof_ids.insert(dof_id); + + n_vars++; + } + } + + opt.set_n_vars(n_vars); + + opt._level_set_sys->solution->close(); +} + + +template +void +MAST::Examples::Bracket3DModel::init_simp_dvs(Opt& opt) { + + libmesh_assert(opt._initialized); + + // + // this assumes that density variable has a constant value per element + // + libmesh_assert_equal_to(opt._density_fetype.family, libMesh::LAGRANGE); + + Real + tol = 1.e-12, + l_frac = 0.4,//_input("length_fraction", "fraction of length along x-axis that is in the bracket", 0.4), + h_frac = 0.4,//_input( "height_fraction", "fraction of length along y-axis that is in the bracket", 0.4), + length = opt._input("length", "length of domain along x-axis", 0.3), + height = opt._input("height", "length of domain along y-axis", 0.3), + x_lim = length * l_frac, + y_lim = height * (1.-h_frac), + frac = opt._input("loadlength_fraction", "fraction of boundary length on which pressure will act", 0.125), + filter_radius = opt._input("filter_radius", "radius of geometric filter for level set field", 0.015); + + unsigned int + sys_num = opt._density_sys->number(), + dof_id = 0, + n_vars = 0; + + Real + val = 0.; + + // + // all ranks will have DVs defined for all variables. So, we should be + // operating on a replicated mesh + // + libmesh_assert(opt._mesh->is_replicated()); + + std::vector local_phi(opt._density_sys->solution->size()); + opt._density_sys->solution->localize(local_phi); + + // iterate over all the element values + libMesh::MeshBase::const_node_iterator + it = opt._mesh->nodes_begin(), + end = opt._mesh->nodes_end(); + + // + // maximum number of dvs is the number of nodes on the level set function + // mesh. We will evaluate the actual number of dvs + // + opt._dv_params.reserve(opt._mesh->n_elem()); + + for ( ; it!=end; it++) { + + const libMesh::Node& n = **it; + + dof_id = n.dof_number(sys_num, 0, 0); + + if ((n(1)-filter_radius) <= y_lim && (n(0)+filter_radius) >= length*(1.-frac)) { + + // + // set value at the constrained points to a small positive number + // material here + // + if (dof_id >= opt._density_sys->solution->first_local_index() && + dof_id < opt._density_sys->solution->last_local_index()) + opt._density_sys->solution->set(dof_id, 1.e0); + } + else { + + std::ostringstream oss; + oss << "dv_" << n_vars; + val = local_phi[dof_id]; + + // + // on the boundary, set everything to be zero, so that there + // is always a boundary there that the optimizer can move + // + if (n(0) < tol || // left boundary + std::fabs(n(0) - length) < tol || // right boundary + std::fabs(n(1) - height) < tol || // top boundary + (n(0) >= x_lim && n(1) <= y_lim)) { + + if (dof_id >= opt._density_sys->solution->first_local_index() && + dof_id < opt._density_sys->solution->last_local_index()) + opt._density_sys->solution->set(dof_id, opt._rho_min); + val = opt._rho_min; + } + + opt._dv_params.push_back(std::pair()); + opt._dv_params[n_vars].first = dof_id; + opt._dv_params[n_vars].second = new MAST::LevelSetParameter(oss.str(), val, &n); + opt._dv_params[n_vars].second->set_as_topology_parameter(true); + opt._dv_dof_ids.insert(dof_id); + + n_vars++; + } + } + + opt.set_n_vars(n_vars); + opt._density_sys->solution->close(); + +} + + + +template +void +MAST::Examples::Bracket3DModel:: +_delete_elems_from_bracket_mesh(Opt& opt, + libMesh::MeshBase &mesh) { + + Real + tol = 1.e-12, + x = -1., + y = -1., + l_frac = 0.4, + w_frac = 0.4, + length = opt._input("length", "length of domain along x-axis", 0.3), + height = opt._input("height", "length of domain along y-axis", 0.3), + x_lim = length * l_frac, + y_lim = height * (1.-w_frac); + + // + // now, remove elements that are outside of the L-bracket domain + // + libMesh::MeshBase::element_iterator + e_it = mesh.elements_begin(), + e_end = mesh.elements_end(); + + for ( ; e_it!=e_end; e_it++) { + + libMesh::Elem* elem = *e_it; + x = length; + y = 0.; + for (unsigned int i=0; in_nodes(); i++) { + const libMesh::Node& n = elem->node_ref(i); + if (x > n(0)) x = n(0); + if (y < n(1)) y = n(1); + } + + // + // delete element if the lowest x,y locations are outside of the bracket + // domain + // + if (x >= x_lim && y<= y_lim) + mesh.delete_elem(elem); + } + + mesh.prepare_for_use(); + + // + // add the two additional boundaries to the boundary info so that + // we can apply loads on them + // + bool + facing_right = false, + facing_down = false; + + e_it = mesh.elements_begin(); + e_end = mesh.elements_end(); + + for ( ; e_it != e_end; e_it++) { + + libMesh::Elem* elem = *e_it; + + if (!elem->on_boundary()) continue; + + for (unsigned int i=0; in_sides(); i++) { + + if (elem->neighbor_ptr(i)) continue; + + std::unique_ptr s(elem->side_ptr(i).release()); + + const libMesh::Point p = s->centroid(); + + facing_right = true; + facing_down = true; + for (unsigned int j=0; jn_nodes(); j++) { + const libMesh::Node& n = s->node_ref(j); + + if (n(0) < x_lim || n(1) > y_lim) { + facing_right = false; + facing_down = false; + } + else if (std::fabs(n(0) - p(0)) > tol) + facing_right = false; + else if (std::fabs(n(1) - p(1)) > tol) + facing_down = false; + } + + if (facing_right) mesh.boundary_info->add_side(elem, i, 6); + if (facing_down) mesh.boundary_info->add_side(elem, i, 7); + } + } + + mesh.boundary_info->sideset_name(6) = "facing_right"; + mesh.boundary_info->sideset_name(7) = "facing_down"; +} + + + +template +void +MAST::Examples::Bracket3DModel::initialize_level_set_solution(Opt& opt) { + + Real + length = opt._input("length", "length of domain along x-axis", 0.3), + height = opt._input("height", "length of domain along y-axis", 0.3), + width = opt._input("width", "length of domain along z-axis", 0.3); + + unsigned int + nx_h = opt._input("initial_level_set_n_holes_in_x", + "number of holes along x-direction for initial level-set field", 6), + ny_h = opt._input("initial_level_set_n_holes_in_y", + "number of holes along y-direction for initial level-set field", 6), + nz_h = opt._input("initial_level_set_n_holes_in_z", + "number of holes along z-direction for initial level-set field", 6), + nx_m = opt._input("level_set_nx_divs", "number of elements of level-set mesh along x-axis", 10), + ny_m = opt._input("level_set_ny_divs", "number of elements of level-set mesh along y-axis", 10), + nz_m = opt._input("level_set_nz_divs", "number of elements of level-set mesh along z-axis", 10); + + MAST::Examples::LevelSetNucleationFunction3D + phi(0., 0., 0., length, height, width, nx_m, ny_m, nz_m, nx_h, ny_h, nz_h); + + opt._level_set_sys_init->initialize_solution(phi); +} + + +#endif // __mast_topology_3d_bracket_model__ diff --git a/examples/structural/base/eyebar_2d_model.h b/examples/structural/base/eyebar_2d_model.h index ddbecc7d..06afaca5 100644 --- a/examples/structural/base/eyebar_2d_model.h +++ b/examples/structural/base/eyebar_2d_model.h @@ -51,11 +51,14 @@ class DisciplineBase; class BoundaryConditionBase; class FunctionBase; class Parameter; +class Solid2DSectionElementPropertyCard; namespace Examples { struct Eyebar2DModel { + using SectionPropertyCardType = MAST::Solid2DSectionElementPropertyCard; + template static Real reference_volume(Opt& opt); @@ -541,7 +544,7 @@ MAST::Examples::Eyebar2DModel::initialize_level_set_solution(Opt& opt) { nx_m = opt._input("level_set_nx_divs", "number of elements of level-set mesh along x-axis", 10), ny_m = opt._input("level_set_ny_divs", "number of elements of level-set mesh along y-axis", 10); - MAST::Examples::LevelSetNucleationFunction + MAST::Examples::LevelSetNucleationFunction2D phi(-0.5*height, -0.5*height, length, height, nx_m, ny_m, nx_h, ny_h); opt._level_set_sys_init->initialize_solution(phi); diff --git a/examples/structural/base/inplane_2d_model.h b/examples/structural/base/inplane_2d_model.h index dc033d66..d0291bb2 100644 --- a/examples/structural/base/inplane_2d_model.h +++ b/examples/structural/base/inplane_2d_model.h @@ -48,12 +48,15 @@ class DisciplineBase; class BoundaryConditionBase; class FunctionBase; class Parameter; +class Solid2DSectionElementPropertyCard; namespace Examples { struct Inplane2DModel { + using SectionPropertyCardType = MAST::Solid2DSectionElementPropertyCard; + template static Real reference_volume(Opt& opt); @@ -497,7 +500,7 @@ MAST::Examples::Inplane2DModel::initialize_level_set_solution(Opt& opt) { nx_m = opt._input("level_set_nx_divs", "number of elements of level-set mesh along x-axis", 10), ny_m = opt._input("level_set_ny_divs", "number of elements of level-set mesh along y-axis", 10); - MAST::Examples::LevelSetNucleationFunction + MAST::Examples::LevelSetNucleationFunction2D phi(0., 0., length, height, nx_m, ny_m, nx_h, ny_h); opt._level_set_sys_init->initialize_solution(phi); diff --git a/examples/structural/base/level_set_nucleation.h b/examples/structural/base/level_set_nucleation.h index 5d305aa6..6d25b735 100644 --- a/examples/structural/base/level_set_nucleation.h +++ b/examples/structural/base/level_set_nucleation.h @@ -28,18 +28,18 @@ namespace MAST { namespace Examples { -class LevelSetNucleationFunction: +class LevelSetNucleationFunction2D: public MAST::FieldFunction { public: - LevelSetNucleationFunction(Real x0, - Real y0, - Real l1, - Real l2, - Real nx_mesh, - Real ny_mesh, - Real nx_holes, - Real ny_holes): + LevelSetNucleationFunction2D(Real x0, + Real y0, + Real l1, + Real l2, + Real nx_mesh, + Real ny_mesh, + Real nx_holes, + Real ny_holes): MAST::FieldFunction("Phi"), _x0 (x0), _y0 (y0), @@ -64,9 +64,9 @@ public MAST::FieldFunction { for (unsigned int i=0; i<_ny_holes; i++) _y_axis_hole_locations.insert(_y0+(i+0.5)*dx); } - - virtual ~LevelSetNucleationFunction() {} - + + virtual ~LevelSetNucleationFunction2D() {} + virtual void operator()(const libMesh::Point& p, const Real t, RealVectorX& v) const { @@ -136,6 +136,145 @@ public MAST::FieldFunction { }; +class LevelSetNucleationFunction3D: +public MAST::FieldFunction { + +public: + LevelSetNucleationFunction3D(Real x0, + Real y0, + Real z0, + Real l1, + Real l2, + Real l3, + Real nx_mesh, + Real ny_mesh, + Real nz_mesh, + Real nx_holes, + Real ny_holes, + Real nz_holes): + MAST::FieldFunction("Phi"), + _x0 (x0), + _y0 (y0), + _z0 (z0), + _l1 (l1), + _l2 (l2), + _l3 (l3), + _nx_mesh (nx_mesh), + _ny_mesh (ny_mesh), + _nz_mesh (nz_mesh), + _nx_holes (nx_holes), + _ny_holes (ny_holes), + _nz_holes (nz_holes), + _pi (acos(-1.)) { + + Real + dx = _l1/(1.*_nx_holes); + + for (unsigned int i=0; i<_nx_holes; i++) + _x_axis_hole_locations.insert(_x0+(i+.5)*dx); + + // + // now, along the y-axis + // + dx = _l2/(1.*_ny_holes); + for (unsigned int i=0; i<_ny_holes; i++) + _y_axis_hole_locations.insert(_y0+(i+0.5)*dx); + + // + // now, along the z-axis + // + dx = _l3/(1.*_nz_holes); + for (unsigned int i=0; i<_nz_holes; i++) + _z_axis_hole_locations.insert(_z0+(i+0.5)*dx); + } + + virtual ~LevelSetNucleationFunction3D() {} + + virtual void operator()(const libMesh::Point& p, + const Real t, + RealVectorX& v) const { + + libmesh_assert_less_equal(t, 1); + libmesh_assert_equal_to(v.size(), 1); + + // + // the libMesh solution projection routine for Lagrange elements + // will query the function value at the nodes. So, we figure + // out which nodes should have zero values set to them. + // if there is one hole in any direction, it will be in the + // center of the domain. If there are more than 1, then two of + // the holes will be on the boundary and others will fill the + // interior evenly. + // + const Real + dx_mesh = _l1/(1.*_nx_holes), + dy_mesh = _l2/(1.*_ny_holes), + dz_mesh = _l3/(1.*_nz_holes); + + std::set::const_iterator + x_it_low = _x_axis_hole_locations.lower_bound(p(0)-dx_mesh), + y_it_low = _y_axis_hole_locations.lower_bound(p(1)-dy_mesh), + z_it_low = _z_axis_hole_locations.lower_bound(p(2)-dz_mesh); + + unsigned int + n = 0; + // + // see if the x-location needs a hole + // + for ( ; x_it_low != _x_axis_hole_locations.end(); x_it_low++) { + if (std::fabs(*x_it_low - p(0)) <= dx_mesh*0.25) { + n++; + break; + } + } + + // + // now check the y-location + // + for ( ; y_it_low != _y_axis_hole_locations.end(); y_it_low++) { + if (std::fabs(*y_it_low - p(1)) <= dy_mesh*0.25) { + n++; + break; + } + } + + // + // now check the z-location + // + for ( ; z_it_low != _z_axis_hole_locations.end(); z_it_low++) { + if (std::fabs(*z_it_low - p(2)) <= dz_mesh*0.25) { + n++; + break; + } + } + + if (n == 3) + v(0) = -1.e0; + else + v(0) = 1.e0; + } + + +protected: + Real + _x0, + _y0, + _z0, + _l1, + _l2, + _l3, + _nx_mesh, + _ny_mesh, + _nz_mesh, + _nx_holes, + _ny_holes, + _nz_holes, + _pi; + std::set _x_axis_hole_locations; + std::set _y_axis_hole_locations; + std::set _z_axis_hole_locations; +}; + } } diff --git a/examples/structural/base/truss_2d_model.h b/examples/structural/base/truss_2d_model.h index 15ea90b6..e4047859 100644 --- a/examples/structural/base/truss_2d_model.h +++ b/examples/structural/base/truss_2d_model.h @@ -49,11 +49,14 @@ class DisciplineBase; class BoundaryConditionBase; class FunctionBase; class Parameter; +class Solid2DSectionElementPropertyCard; namespace Examples { struct Truss2DModel { + using SectionPropertyCardType = MAST::Solid2DSectionElementPropertyCard; + template static Real reference_volume(Opt& opt); @@ -457,7 +460,7 @@ MAST::Examples::Truss2DModel::initialize_level_set_solution(Opt& opt) { nx_m = opt._input("level_set_nx_divs", "number of elements of level-set mesh along x-axis", 10), ny_m = opt._input("level_set_ny_divs", "number of elements of level-set mesh along y-axis", 10); - MAST::Examples::LevelSetNucleationFunction + MAST::Examples::LevelSetNucleationFunction2D phi(0., 0., length, height, nx_m, ny_m, nx_h, ny_h); opt._level_set_sys_init->initialize_solution(phi); diff --git a/examples/structural/example_5/example_5.cpp b/examples/structural/example_5/example_5.cpp index 736262f2..cbbad772 100644 --- a/examples/structural/example_5/example_5.cpp +++ b/examples/structural/example_5/example_5.cpp @@ -202,8 +202,8 @@ public MAST::FunctionEvaluation { MAST::FilterBase* _filter; - MAST::MaterialPropertyCardBase *_m_card1, *_m_card2; - MAST::ElementPropertyCardBase *_p_card1, *_p_card2; + MAST::MaterialPropertyCardBase *_m_card; + MAST::ElementPropertyCardBase *_p_card; PhiMeshFunction* _level_set_function; MAST::LevelSetBoundaryVelocity* _level_set_vel; @@ -369,20 +369,12 @@ public MAST::FunctionEvaluation { _field_functions.insert(k_f); _field_functions.insert(cp_f); - _m_card1 = new MAST::IsotropicMaterialPropertyCard; - _m_card2 = new MAST::IsotropicMaterialPropertyCard; - _m_card1->add(*E_f); - _m_card1->add(*rho_f); - _m_card1->add(*nu_f); - _m_card1->add(*k_f); - _m_card1->add(*cp_f); - - // material for void - _m_card2->add(*E_v_f); - _m_card2->add(*rho_f); - _m_card2->add(*nu_f); - _m_card2->add(*k_f); - _m_card2->add(*cp_f); + _m_card = new MAST::IsotropicMaterialPropertyCard; + _m_card->add(*E_f); + _m_card->add(*rho_f); + _m_card->add(*nu_f); + _m_card->add(*k_f); + _m_card->add(*cp_f); } @@ -416,44 +408,26 @@ public MAST::FunctionEvaluation { _field_functions.insert(kappa_f); _field_functions.insert(hoff_f); - MAST::Solid2DSectionElementPropertyCard - *p_card1 = new MAST::Solid2DSectionElementPropertyCard, - *p_card2 = new MAST::Solid2DSectionElementPropertyCard; - - _p_card1 = p_card1; - _p_card2 = p_card2; + typename T::SectionPropertyCardType + *p_card = new typename T::SectionPropertyCardType; + + _p_card = p_card; // set nonlinear strain if requested bool nonlinear = _input("if_nonlinear", "flag to turn on/off nonlinear strain", false); - if (nonlinear) _p_card1->set_strain(MAST::NONLINEAR_STRAIN); - _p_card2->set_strain(MAST::LINEAR_STRAIN); - - p_card1->add(*th_f); - p_card1->add(*kappa_f); - p_card1->add(*hoff_f); - p_card1->set_material(*_m_card1); + if (nonlinear) _p_card->set_strain(MAST::NONLINEAR_STRAIN); // property card for void - p_card2->add(*th_f); - p_card2->add(*kappa_f); - p_card2->add(*hoff_f); - p_card2->set_material(*_m_card2); - - _discipline->set_property_for_subdomain(0, *p_card1); - _discipline->set_property_for_subdomain(1, *p_card1); + p_card->add(*th_f); + p_card->add(*kappa_f); + p_card->add(*hoff_f); + p_card->set_material(*_m_card); - // inactive - _discipline->set_property_for_subdomain(3, *p_card2); + _discipline->set_property_for_subdomain(0, *p_card); - // negative level set - _discipline->set_property_for_subdomain(6, *p_card2); - _discipline->set_property_for_subdomain(7, *p_card2); - _indicator_discipline->set_property_for_subdomain(0, *p_card1); - _indicator_discipline->set_property_for_subdomain(1, *p_card1); - _indicator_discipline->set_property_for_subdomain(3, *p_card2); - _indicator_discipline->set_property_for_subdomain(4, *p_card2); + _indicator_discipline->set_property_for_subdomain(0, *p_card); } // @@ -1455,10 +1429,8 @@ public MAST::FunctionEvaluation { _indicator_discipline (nullptr), _level_set_discipline (nullptr), _filter (nullptr), - _m_card1 (nullptr), - _m_card2 (nullptr), - _p_card1 (nullptr), - _p_card2 (nullptr), + _m_card (nullptr), + _p_card (nullptr), _level_set_function (nullptr), _level_set_vel (nullptr), _output (nullptr), @@ -1494,12 +1466,6 @@ public MAST::FunctionEvaluation { _nsp = new MAST::StructuralNearNullVectorSpace; _sys->nonlinear_solver->nearnullspace_object = _nsp; - // - // ask structure to use Mindlin bending operator - // - dynamic_cast(*_p_card1).set_bending_model(MAST::MINDLIN); - dynamic_cast(*_p_card2).set_bending_model(MAST::MINDLIN); - ///////////////////////////////////////////////// // now initialize the design data. ///////////////////////////////////////////////// @@ -1596,10 +1562,8 @@ public MAST::FunctionEvaluation { delete _nsp; - delete _m_card1; - delete _m_card2; - delete _p_card1; - delete _p_card2; + delete _m_card; + delete _p_card; delete _eq_sys; delete _mesh_refinement; diff --git a/examples/structural/example_6/example_6.cpp b/examples/structural/example_6/example_6.cpp index 3d885464..1f21ff32 100644 --- a/examples/structural/example_6/example_6.cpp +++ b/examples/structural/example_6/example_6.cpp @@ -44,10 +44,12 @@ #include "solver/slepc_eigen_solver.h" #include "property_cards/isotropic_material_property_card.h" #include "property_cards/solid_2d_section_element_property_card.h" +#include "property_cards/isotropic_element_property_card_3D.h" #include "optimization/gcmma_optimization_interface.h" #include "optimization/npsol_optimization_interface.h" #include "optimization/function_evaluation.h" #include "examples/structural/base/bracket_2d_model.h" +#include "examples/structural/base/bracket_3d_model.h" #include "examples/structural/base/inplane_2d_model.h" #include "examples/structural/base/truss_2d_model.h" #include "examples/structural/base/eyebar_2d_model.h" @@ -361,8 +363,8 @@ public MAST::FunctionEvaluation { _field_functions.insert(kappa_f); _field_functions.insert(hoff_f); - MAST::Solid2DSectionElementPropertyCard - *p_card = new MAST::Solid2DSectionElementPropertyCard; + typename T::SectionPropertyCardType + *p_card = new typename T::SectionPropertyCardType; _p_card = p_card; @@ -463,14 +465,18 @@ public MAST::FunctionEvaluation { // this will create a localized vector in _level_set_sys->curret_local_solution _density_sys->update(); - // create a serialized vector for use in interpolation + // create a localized vector for use in interpolation std::unique_ptr> - serial_density_sol(libMesh::NumericVector::build(_sys->comm()).release()); - serial_density_sol->init(_density_sys->solution->size(), false, libMesh::SERIAL); - _density_sys->solution->localize(*serial_density_sol); + local_density_sol(libMesh::NumericVector::build(_sys->comm()).release()); + local_density_sol->init(_density_sys->n_dofs(), + _density_sys->n_local_dofs(), + _density_sys->get_dof_map().get_send_list(), + false, + libMesh::GHOSTED); + _density_sys->solution->localize(*local_density_sol); _density_function->clear(); - _density_function->init(*serial_density_sol, true); + _density_function->init(*local_density_sol, true); _sys->solution->zero(); ////////////////////////////////////////////////////////////////////// @@ -552,7 +558,7 @@ public MAST::FunctionEvaluation { ////////////////////////////////////////////////////////////////////// // evaluate the volume for used in the problem setup - _evaluate_volume(*serial_density_sol, sens_vecs, &vol, nullptr); + _evaluate_volume(*local_density_sol, sens_vecs, &vol, nullptr); libMesh::out << "volume: " << vol << std::endl; // evaluate the output based on specified problem type @@ -599,7 +605,7 @@ public MAST::FunctionEvaluation { } else if (_problem == "volume_stress") { - _evaluate_volume(*serial_density_sol, sens_vecs, nullptr, &obj_grad); + _evaluate_volume(*local_density_sol, sens_vecs, nullptr, &obj_grad); for (unsigned int i=0; i::build(_sys->comm()).release(); - vec->init(_density_sys->solution->size(), false, libMesh::SERIAL); + vec->init(_density_sys->n_dofs(), + _density_sys->n_local_dofs(), + _density_sys->get_dof_map().get_send_list(), + false, + libMesh::GHOSTED); _filter->compute_filtered_values(nonzero_val, *vec, false); dphi_vecs[i] = vec; @@ -674,7 +684,7 @@ public MAST::FunctionEvaluation { for ( unsigned int i=0; i<_n_vars; i++) dphi_vecs[i]->close(); - // we will use this serialized vector to initialize the mesh function, + // we will use this ghosted vector to initialize the mesh function, // which is setup to reuse this vector, so we have to store it for ( unsigned int i=0; i<_n_vars; i++) _density_function->init_sens(*_dv_params[i].second, *dphi_vecs[i], true); @@ -957,13 +967,17 @@ public MAST::FunctionEvaluation { base_phi.set(_dv_params[i].first, x[i]); base_phi.close(); _filter->compute_filtered_values(base_phi, *_density_sys->solution); - // create a serialized vector for use in interpolation + // create a ghosted vector for use in interpolation std::unique_ptr> - serial_density_sol(libMesh::NumericVector::build(_sys->comm()).release()); - serial_density_sol->init(_density_sys->solution->size(), false, libMesh::SERIAL); - _density_sys->solution->localize(*serial_density_sol); + local_density_sol(libMesh::NumericVector::build(_sys->comm()).release()); + local_density_sol->init(_density_sys->n_dofs(), + _density_sys->n_local_dofs(), + _density_sys->get_dof_map().get_send_list(), + false, + libMesh::GHOSTED); + _density_sys->solution->localize(*local_density_sol); - _density_function->init(*serial_density_sol, true); + _density_function->init(*local_density_sol, true); std::vector eval_grads(this->n_ineq(), false); std::vector f(this->n_ineq(), 0.), grads; @@ -1079,18 +1093,13 @@ public MAST::FunctionEvaluation { // density function is used by elasticity modulus function. So, we // initialize this here - _density_function = new MAST::MeshFieldFunction(*_density_sys, "rho", libMesh::SERIAL); + _density_function = new MAST::MeshFieldFunction(*_density_sys, "rho", libMesh::GHOSTED); _init_material(); T::init_structural_loads(*this); _init_section_property(); _initialized = true; - // - // ask structure to use Mindlin bending operator - // - dynamic_cast(*_p_card).set_bending_model(MAST::MINDLIN); - ///////////////////////////////////////////////// // now initialize the design data. ///////////////////////////////////////////////// @@ -1390,6 +1399,11 @@ int main(int argc, char* argv[]) { (new TopologyOptimizationSIMP (init.comm(), input)); } + else if (mesh == "bracket3d") { + top_opt.reset + (new TopologyOptimizationSIMP + (init.comm(), input)); + } else libmesh_error(); diff --git a/examples/structural/example_8/example_8.cpp b/examples/structural/example_8/example_8.cpp index 8122c28c..6f5183cd 100644 --- a/examples/structural/example_8/example_8.cpp +++ b/examples/structural/example_8/example_8.cpp @@ -249,8 +249,8 @@ public MAST::FunctionEvaluation { MAST::FilterBase* _filter; - MAST::MaterialPropertyCardBase *_m_card1, *_m_card2; - MAST::ElementPropertyCardBase *_p_card1, *_p_card2; + MAST::MaterialPropertyCardBase *_m_card; + MAST::ElementPropertyCardBase *_p_card; PhiMeshFunction* _level_set_function; libMesh::ExodusII_IO* _output; @@ -399,20 +399,12 @@ public MAST::FunctionEvaluation { _field_functions.insert(k_f); _field_functions.insert(cp_f); - _m_card1 = new MAST::IsotropicMaterialPropertyCard; - _m_card2 = new MAST::IsotropicMaterialPropertyCard; - _m_card1->add(*_Ef); - _m_card1->add(*rho_f); - _m_card1->add(*nu_f); - _m_card1->add(*k_f); - _m_card1->add(*cp_f); - - // material for void - _m_card2->add(*_Ef); - _m_card2->add(*rho_f); - _m_card2->add(*nu_f); - _m_card2->add(*k_f); - _m_card2->add(*cp_f); + _m_card = new MAST::IsotropicMaterialPropertyCard; + _m_card->add(*_Ef); + _m_card->add(*rho_f); + _m_card->add(*nu_f); + _m_card->add(*k_f); + _m_card->add(*cp_f); } @@ -446,39 +438,22 @@ public MAST::FunctionEvaluation { _field_functions.insert(kappa_f); _field_functions.insert(hoff_f); - MAST::Solid2DSectionElementPropertyCard - *p_card1 = new MAST::Solid2DSectionElementPropertyCard, - *p_card2 = new MAST::Solid2DSectionElementPropertyCard; - - _p_card1 = p_card1; - _p_card2 = p_card2; + typename T::SectionPropertyCardType + *p_card = new typename T::SectionPropertyCardType; + + _p_card = p_card; // set nonlinear strain if requested bool nonlinear = _input("if_nonlinear", "flag to turn on/off nonlinear strain", false); - if (nonlinear) _p_card1->set_strain(MAST::NONLINEAR_STRAIN); - _p_card2->set_strain(MAST::LINEAR_STRAIN); + if (nonlinear) _p_card->set_strain(MAST::NONLINEAR_STRAIN); - p_card1->add(*th_f); - p_card1->add(*hoff_f); - p_card1->add(*kappa_f); - p_card1->set_material(*_m_card1); - - // property card for void - p_card2->add(*th_f); - p_card2->add(*hoff_f); - p_card2->add(*kappa_f); - p_card2->set_material(*_m_card2); - - _discipline->set_property_for_subdomain(0, *p_card1); - _discipline->set_property_for_subdomain(1, *p_card1); - - // inactive - _discipline->set_property_for_subdomain(3, *p_card2); + p_card->add(*th_f); + p_card->add(*hoff_f); + p_card->add(*kappa_f); + p_card->set_material(*_m_card); - // negative level set - _discipline->set_property_for_subdomain(6, *p_card2); - _discipline->set_property_for_subdomain(7, *p_card2); + _discipline->set_property_for_subdomain(0, *p_card); } // @@ -1312,10 +1287,8 @@ public MAST::FunctionEvaluation { _discipline (nullptr), _level_set_discipline (nullptr), _filter (nullptr), - _m_card1 (nullptr), - _m_card2 (nullptr), - _p_card1 (nullptr), - _p_card2 (nullptr), + _m_card (nullptr), + _p_card (nullptr), _level_set_function (nullptr), _output (nullptr) { @@ -1349,12 +1322,6 @@ public MAST::FunctionEvaluation { _nsp = new MAST::StructuralNearNullVectorSpace; _sys->nonlinear_solver->nearnullspace_object = _nsp; - // - // ask structure to use Mindlin bending operator - // - dynamic_cast(*_p_card1).set_bending_model(MAST::MINDLIN); - dynamic_cast(*_p_card2).set_bending_model(MAST::MINDLIN); - ///////////////////////////////////////////////// // now initialize the design data. ///////////////////////////////////////////////// @@ -1435,10 +1402,8 @@ public MAST::FunctionEvaluation { delete _nsp; - delete _m_card1; - delete _m_card2; - delete _p_card1; - delete _p_card2; + delete _m_card; + delete _p_card; delete _eq_sys; delete _mesh_refinement; diff --git a/src/base/assembly_base.cpp b/src/base/assembly_base.cpp index afcf274c..6340dab7 100644 --- a/src/base/assembly_base.cpp +++ b/src/base/assembly_base.cpp @@ -593,14 +593,14 @@ calculate_output_adjoint_sensitivity(const libMesh::NumericVector& X, void MAST::AssemblyBase::calculate_output_adjoint_sensitivity_multiple_parameters_no_direct - (const libMesh::NumericVector& X, - bool if_localize_sol, - const libMesh::NumericVector& dq_dX, - const std::vector& p_vec, - MAST::AssemblyElemOperations& elem_ops, - MAST::OutputAssemblyElemOperations& output, - std::vector& sens) { - +(const libMesh::NumericVector& X, + bool if_localize_sol, + const libMesh::NumericVector& dq_dX, + const std::vector& p_vec, + MAST::AssemblyElemOperations& elem_ops, + MAST::OutputAssemblyElemOperations& output, + std::vector& sens) { + libmesh_assert(_discipline); libmesh_assert(_system); libmesh_assert_equal_to(sens.size(), p_vec.size()); @@ -613,7 +613,7 @@ MAST::AssemblyBase::calculate_output_adjoint_sensitivity_multiple_parameters_no_ // add vectors before computing sensitivity for (unsigned int i=0; i + fe(_elem.init_fe(true, false)); + + const std::vector& JxW = fe->get_JxW(); + const std::vector& xyz = fe->get_xyz(); + const unsigned int + n_phi = (unsigned int)fe->n_shape_functions(), + n1 =6, + n2 =3*n_phi, + n3 =30, + n_nodes =_elem.get_reference_elem().n_nodes(); + + RealMatrixX + material_mat, + mat_x = RealMatrixX::Zero(6,3), + mat_y = RealMatrixX::Zero(6,3), + mat_z = RealMatrixX::Zero(6,3), + mat1_n1n2 = RealMatrixX::Zero(n1, n2), + mat2_n2n2 = RealMatrixX::Zero(n2, n2), + mat3_3n2 = RealMatrixX::Zero(3, n2), + mat4_33 = RealMatrixX::Zero(3, 3), + mat5_n1n3 = RealMatrixX::Zero(n1, n3), + mat6_n2n3 = RealMatrixX::Zero(n2, n3), + mat7_3n3 = RealMatrixX::Zero(3, n3), + Gmat = RealMatrixX::Zero(6, n3), + K_alphaalpha = RealMatrixX::Zero(n3, n3), + K_ualpha = RealMatrixX::Zero(n2, n3), + K_corr = RealMatrixX::Zero(n2, n2); + RealVectorX + strain = RealVectorX::Zero(6), + stress = RealVectorX::Zero(6), + vec1_n1 = RealVectorX::Zero(n1), + vec2_n2 = RealVectorX::Zero(n2), + vec3_3 = RealVectorX::Zero(3), + local_disp= RealVectorX::Zero(n2), + f_alpha = RealVectorX::Zero(n3), + alpha = RealVectorX::Zero(n3);//*_incompatible_sol; + + // copy the values from the global to the local element + local_disp.topRows(n2) = _local_sol.topRows(n2); + + std::unique_ptr > mat_stiff = + _property.stiffness_A_matrix(*this); + + MAST::FEMOperatorMatrix + Bmat_lin, + Bmat_nl_x, + Bmat_nl_y, + Bmat_nl_z, + Bmat_nl_u, + Bmat_nl_v, + Bmat_nl_w, + Bmat_inc; + // six stress components, related to three displacements + Bmat_lin.reinit(n1, 3, n_nodes); + Bmat_nl_x.reinit(3, 3, n_nodes); + Bmat_nl_y.reinit(3, 3, n_nodes); + Bmat_nl_z.reinit(3, 3, n_nodes); + Bmat_nl_u.reinit(3, 3, n_nodes); + Bmat_nl_v.reinit(3, 3, n_nodes); + Bmat_nl_w.reinit(3, 3, n_nodes); + Bmat_inc.reinit(n1, n3, 1); // six stress-strain components + + /////////////////////////////////////////////////////////////////////// + // second for loop to calculate the residual and stiffness contributions + for (unsigned int qp=0; qpderivative(p, xyz[qp], _time, material_mat); + + this->initialize_green_lagrange_strain_operator(qp, + *fe, + local_disp, + strain, + mat_x, mat_y, mat_z, + Bmat_lin, + Bmat_nl_x, + Bmat_nl_y, + Bmat_nl_z, + Bmat_nl_u, + Bmat_nl_v, + Bmat_nl_w); + //this->initialize_incompatible_strain_operator(qp, *_fe, Bmat_inc, Gmat); + + // calculate the stress + stress = material_mat * (strain + Gmat * alpha); + + // residual from incompatible modes + f_alpha += JxW[qp] * Gmat.transpose() * stress; + + // calculate contribution to the residual + // linear strain operator + Bmat_lin.vector_mult_transpose(vec2_n2, stress); + f.topRows(n2) += JxW[qp] * vec2_n2; + + if (_property.strain_type() == MAST::NONLINEAR_STRAIN) { + + // nonlinear strain operator + // x + vec3_3 = mat_x.transpose() * stress; + Bmat_nl_x.vector_mult_transpose(vec2_n2, vec3_3); + f.topRows(n2) += JxW[qp] * vec2_n2; + + // y + vec3_3 = mat_y.transpose() * stress; + Bmat_nl_y.vector_mult_transpose(vec2_n2, vec3_3); + f.topRows(n2) += JxW[qp] * vec2_n2; + + // z + vec3_3 = mat_z.transpose() * stress; + Bmat_nl_z.vector_mult_transpose(vec2_n2, vec3_3); + f.topRows(n2) += JxW[qp] * vec2_n2; + } + + if (request_jacobian) { + + // the strain includes the following expansion + // delta_epsilon = B_lin + mat_x B_x + mat_y B_y + mat_z B_z + // Hence, the tangent stiffness matrix will include + // components from epsilon^T C epsilon + + //////////////////////////////////////////////////////// + // B_lin^T C B_lin + Bmat_lin.left_multiply(mat1_n1n2, material_mat); + Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + if (_property.strain_type() == MAST::NONLINEAR_STRAIN) { + + // B_x^T mat_x^T C B_lin + mat3_3n2 = mat_x.transpose() * mat1_n1n2; + Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + // B_y^T mat_y^T C B_lin + mat3_3n2 = mat_y.transpose() * mat1_n1n2; + Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + // B_z^T mat_z^T C B_lin + mat3_3n2 = mat_z.transpose() * mat1_n1n2; + Bmat_nl_z.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + /////////////////////////////////////////////////////// + for (unsigned int i_dim=0; i_dim<3; i_dim++) { + switch (i_dim) { + case 0: + Bmat_nl_x.left_multiply(mat1_n1n2, mat_x); + break; + + case 1: + Bmat_nl_y.left_multiply(mat1_n1n2, mat_y); + break; + + case 2: + Bmat_nl_z.left_multiply(mat1_n1n2, mat_z); + break; + } + + // B_lin^T C mat_x_i B_x_i + mat1_n1n2 = material_mat * mat1_n1n2; + Bmat_lin.right_multiply_transpose(mat2_n2n2, mat1_n1n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + // B_x^T mat_x^T C mat_x B_x + mat3_3n2 = mat_x.transpose() * mat1_n1n2; + Bmat_nl_x.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + // B_y^T mat_y^T C mat_x B_x + mat3_3n2 = mat_y.transpose() * mat1_n1n2; + Bmat_nl_y.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + // B_z^T mat_z^T C mat_x B_x + mat3_3n2 = mat_z.transpose() * mat1_n1n2; + Bmat_nl_z.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + } + + // use the stress to calculate the final contribution + // to the Jacobian stiffness matrix + mat4_33(0,0) = stress(0); + mat4_33(1,1) = stress(1); + mat4_33(2,2) = stress(2); + mat4_33(0,1) = mat4_33(1,0) = stress(3); + mat4_33(1,2) = mat4_33(2,1) = stress(4); + mat4_33(0,2) = mat4_33(2,0) = stress(5); + + // u-disp + Bmat_nl_u.left_multiply(mat3_3n2, mat4_33); + Bmat_nl_u.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + // v-disp + Bmat_nl_v.left_multiply(mat3_3n2, mat4_33); + Bmat_nl_v.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + + // w-disp + Bmat_nl_w.left_multiply(mat3_3n2, mat4_33); + Bmat_nl_w.right_multiply_transpose(mat2_n2n2, mat3_3n2); + jac.topLeftCorner(n2, n2) += JxW[qp] * mat2_n2n2; + } + } + } + + // if jacobian is requested, add a small diagonal value for the + // rotational dofs + if (request_jacobian) + jac.bottomRightCorner(n2, n2) += RealMatrixX::Identity(n2, n2) * + 1.0e-20 * jac.diagonal().maxCoeff(); + return request_jacobian; } diff --git a/src/level_set/filter_base.cpp b/src/level_set/filter_base.cpp index 718604ef..23629c15 100644 --- a/src/level_set/filter_base.cpp +++ b/src/level_set/filter_base.cpp @@ -102,7 +102,6 @@ MAST::FilterBase::compute_filtered_values(std::map& nonzero_ bool close_vector) const { libmesh_assert_equal_to(output.size(), _filter_map.size()); - libmesh_assert_equal_to(output.type(), libMesh::SERIAL); output.zero(); @@ -119,10 +118,15 @@ MAST::FilterBase::compute_filtered_values(std::map& nonzero_ for ( ; vec_it != vec_end; vec_it++) { if (nonzero_input.count(vec_it->first)) { - if (_dv_dof_ids.count(map_it->first)) - output.add(map_it->first, nonzero_input[vec_it->first] * vec_it->second); - else - output.set(map_it->first, nonzero_input[map_it->first]); + if (output.type() == libMesh::SERIAL || + (map_it->first >= output.first_local_index() && + map_it->first < output.last_local_index())) { + + if (_dv_dof_ids.count(map_it->first)) + output.add(map_it->first, nonzero_input[vec_it->first] * vec_it->second); + else + output.set(map_it->first, nonzero_input[map_it->first]); + } } } } diff --git a/src/level_set/filter_base.h b/src/level_set/filter_base.h index ca360d69..e918f8fc 100644 --- a/src/level_set/filter_base.h +++ b/src/level_set/filter_base.h @@ -56,8 +56,9 @@ namespace MAST { /*! * for large problems it is more efficient to specify only the non-zero entries in the input vector in - * \p nonzero_vals. Here, \p output is expected to be of type SERIAL vector. All ranks in the + * \p nonzero_vals. If \p output is of type SERIAL then all ranks in the * communicator will perform the same operaitons and provide an identical \p output vector. + * If \p output is of type PARALLEL or GHOSTED then only the local contributions will be set. * If \p close_vector is \p true then \p output.close() will be called in this * routines, otherwise not. */