Feature/add beam library

examples/structural/CMakeLists.txt

@@ -6,3 +6,4 @@ add_subdirectory(example_5) # topology optimization
 add_subdirectory(example_6) # SIMP topology optimization
 add_subdirectory(example_7) # NastranIO input for Nastran BDF mesh
 add_subdirectory(example_8) # Homogenized level-set topology optimization
+add_subdirectory(example_9) # Cross section property analysis

examples/structural/example_9/CMakeLists.txt

@@ -0,0 +1,22 @@
+add_executable(structural_example_9
+    example_9.cpp)
+
+target_link_libraries(structural_example_9 mast)
+
+install(TARGETS structural_example_9
+    RUNTIME DESTINATION ${CMAKE_INSTALL_PREFIX}/examples)
+
+# Test on single processor, PETSc built-in LU direct linear solver (sequential).
+add_test(NAME structural_example_9
+    COMMAND $<TARGET_FILE:structural_example_9> -ksp_type preonly -pc_type lu -options_view)
+set_tests_properties(structural_example_9
+    PROPERTIES
+        LABELS "SEQ")
+
+# Test multiple processors, parallel direct solver using external MUMPS package.
+add_test(NAME structural_example_9_mpi
+    COMMAND ${MPIEXEC_EXECUTABLE} -np 2 $<TARGET_FILE:structural_example_9>
+            -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package mumps -options_view)
+set_tests_properties(structural_example_9_mpi
+    PROPERTIES
+        LABELS "MPI")

examples/structural/example_9/example_9.cpp

@@ -0,0 +1,194 @@
+// C/C++ Includes
+#include <iostream>
+
+// libMesh Includes
+#include "libmesh/libmesh.h"
+#include "libmesh/replicated_mesh.h"
+#include "libmesh/mesh_generation.h"
+#include "libmesh/point.h"
+#include "libmesh/edge_edge2.h"
+#include "libmesh/equation_systems.h"
+#include "libmesh/exodusII_io.h"
+
+// MAST Includes
+#include "base/nonlinear_system.h"
+#include "elasticity/structural_system_initialization.h"
+#include "base/physics_discipline_base.h"
+#include "boundary_condition/dirichlet_boundary_condition.h"
+#include "base/parameter.h"
+#include "base/constant_field_function.h"
+#include "property_cards/isotropic_material_property_card.h"
+#include "property_cards/solid_1d_section_element_property_card.h"
+#include "base/nonlinear_implicit_assembly.h"
+#include "elasticity/structural_nonlinear_assembly.h"
+
+#include "elasticity/structural_buckling_eigenproblem_assembly.h"
+#include "elasticity/structural_buckling_eigenproblem_elem_operations.h"
+#include "solver/slepc_eigen_solver.h"
+#include <boundary_condition/point_load_condition.h>
+
+#include "property_cards/solid_1d_bar_section_element_property_card.h"
+
+#include "property_cards/cross_section_property_pilkey.h"
+
+#include "libfort/fort.hpp"
+
+
+// This is the main function which runs when the compiled code is called.
+int main(int argc, const char** argv)
+{
+    // Initialize llibMesh Library
+    libMesh::LibMeshInit init(argc, argv);
+
+    // The target number of elements to be used in the mesh of the cross 
+    // section.  The section is meshed with triangle which only allows you to
+    // specified a desired area for an individual element and not a number of 
+    // elements, so the number of elements in the section mesh won't match this
+    // exactly.
+    const uint n_target_elements = 3000;
+
+    // Define Material Properties as MAST Parameters
+    // For a cross section property analysis, only Poisson's ratio is required,
+    // although this could change if a cross section has different properties
+    // throughout the section (i.e. a composite section).
+    MAST::Parameter nu("nu_param", 0.33);             // Poisson's ratio
+
+    // Define Section Properties as MAST Parameters
+    MAST::Parameter DIM1("DIM1", 3.234);
+    MAST::Parameter DIM2("DIM2", 1.422);
+    MAST::Parameter offset_y("offy_param", 0.587);  // Section offset in y-direction
+    MAST::Parameter offset_z("offz_param", -1.054); // Section offset in z-direction
+
+    // Define Sensitivity Parameters
+    std::vector<MAST::Parameter*> sens_params = {&DIM1, &DIM2};
+    uint n_s = sens_params.size();
+
+    // Create field functions to dsitribute these constant parameters throughout the model
+    // Section Property Field Functions
+    MAST::ConstantFieldFunction DIM1_f("DIM1", DIM1);
+    MAST::ConstantFieldFunction DIM2_f("DIM2", DIM2);
+    MAST::ConstantFieldFunction offsety_f("hy_off", offset_y);
+    MAST::ConstantFieldFunction offsetz_f("hz_off", offset_z);
+    // Material Property Field Functions
+    MAST::ConstantFieldFunction nu_f("nu", nu);
+
+    // Initialize the material
+    MAST::IsotropicMaterialPropertyCard material;                   
+
+    // Add the material property constant field functions to the material card
+    material.add(nu_f);
+
+    // Initialize the section
+    MAST::Solid1DBarSectionElementPropertyCard section;
+
+    // Add the section property constant field functions to the section card
+    section.add(DIM1_f);
+    section.add(DIM2_f);
+    section.add(offsety_f);
+    section.add(offsetz_f);
+
+    // Add the material card to the section card
+    section.set_material(material);
+
+    // Specify a section orientation point and add it to the section.
+    RealVectorX orientation = RealVectorX::Zero(3);
+    orientation(1) = 1.0;
+    section.y_vector() = orientation;
+
+    // Now initialize the section
+    section.init(init, n_target_elements);
+
+    const libMesh::Point point(4.3, -3.5, -6.7);
+    const Real time = 8.22;
+
+    // Setup the table for section property console output
+    fort::table properties_out;
+    properties_out << fort::header << "Name" << "Property" << "Value" << fort::endr;
+
+    fort::table dproperties_out;
+    dproperties_out << fort::header << "Property" << "Value" << fort::endr;
+
+    // Get Area
+    const MAST::FieldFunction<Real>& Area = section.A();
+    Real A, dA;
+    Area(point, time, A);
+    Area.derivative(DIM1, point, time, dA);
+    properties_out << "Area" << "A" << std::to_string(A) << fort::endr;
+    dproperties_out << "dA" << std::to_string(dA) << fort::endr;
+
+    // Get Second Area Moments
+    const MAST::FieldFunction<RealMatrixX>& SecondAreaMoments = section.I();
+    Real Izz, Iyy, Izy;
+    RealMatrixX I, dI;
+    SecondAreaMoments(point, time, I);
+    SecondAreaMoments.derivative(DIM1, point, time, dI);
+    Izz = I(0,0);
+    Iyy = I(1,1);
+    Izy = I(0,1);
+    properties_out << "Inertia" << "I_zz" << std::to_string(Izz) << fort::endr;
+    properties_out << "Inertia" << "I_yy" << std::to_string(Iyy) << fort::endr;
+    properties_out << "Inertia" << "I_zy" << std::to_string(Izy) << fort::endr;
+    dproperties_out << "dI_zz" << std::to_string(dI(0,0)) << fort::endr;
+    dproperties_out << "dI_yy" << std::to_string(dI(1,1)) << fort::endr;
+    dproperties_out << "dI_zy" << std::to_string(dI(0,1)) << fort::endr;
+
+    // Get Second Area Polar Moment
+    const MAST::FieldFunction<Real>& PolarInteria = section.Ip();
+    Real Ip, dIp;
+    PolarInteria(point, time, Ip);
+    PolarInteria.derivative(DIM1, point, time, dIp);
+    properties_out << "Polar Inertia" << "I_xx" << std::to_string(Ip) << fort::endr;
+    dproperties_out << "dI_xx" << std::to_string(dIp) << fort::endr;
+
+    // Get Torsion Constant
+    MAST::FieldFunction<Real>& TorsionConstant = section.J();
+    Real J, dJ;
+    TorsionConstant(point, time, J);
+    TorsionConstant.derivative(DIM1, point, time, dJ);
+    properties_out << "Torsion Constant" << "J" << std::to_string(J) << fort::endr;
+    dproperties_out << "dJ" << std::to_string(dJ) << fort::endr;
+
+    // Get first area moments
+    const MAST::FieldFunction<Real>& AreaMomentY = section.Ay();
+    const MAST::FieldFunction<Real>& AreaMomentZ = section.Az();
+    Real Az, Ay, dAy, dAz;
+    AreaMomentY(point, time, Ay);
+    AreaMomentY.derivative(DIM1, point, time, dAy);
+    AreaMomentZ(point, time, Az);
+    AreaMomentZ.derivative(DIM1, point, time, dAz);
+    properties_out << "First Area Moment" << "A_z" << std::to_string(Ay) << fort::endr;
+    dproperties_out << "dA_z" << std::to_string(dAz) << fort::endr;
+    properties_out << "First Area Moment" << "A_y" << std::to_string(Az) << fort::endr;
+    dproperties_out << "dA_y" << std::to_string(dAy) << fort::endr;
+
+    // Shear Coefficeints
+    MAST::FieldFunction<RealMatrixX>& ShearCoefficients = section.Kap();
+    RealMatrixX Kappa, dKappa;
+    ShearCoefficients(point, time, Kappa);
+    ShearCoefficients.derivative(DIM1, point, time, dKappa);
+    properties_out << "Shear Coefficient" << "kappa_zz" << std::to_string(Kappa(0,0)) << fort::endr;
+    properties_out << "Shear Coefficient" << "kappa_yy" << std::to_string(Kappa(1,1)) << fort::endr;
+    properties_out << "Shear Coefficient" << "kappa_zy" << std::to_string(Kappa(0,1)) << fort::endr;
+    dproperties_out << "dkappa_zz" << std::to_string(dKappa(0,0)) << fort::endr;
+    dproperties_out << "dkappa_yy" << std::to_string(dKappa(1,1)) << fort::endr;
+    dproperties_out << "dkappa_zy" << std::to_string(dKappa(0,1)) << fort::endr;
+
+    // Get warping constant
+    MAST::FieldFunction<Real>& WarpingConstant = section.Gam();
+    Real W, dW;
+    WarpingConstant(point, time, W);
+    WarpingConstant.derivative(DIM1, point, time, dW);
+    properties_out << "Warping Constant" << "W" << std::to_string(W) << fort::endr;
+    dproperties_out << "dW" << std::to_string(dW) << fort::endr;
+
+    // Get shear center
+    const libMesh::Point shear_center = section.get_shear_center(point, time);
+    properties_out << "Shear Center" << "z_sc" << std::to_string(shear_center(0)) << fort::endr;
+    properties_out << "Shear Center" << "y_sc" << std::to_string(shear_center(1)) << fort::endr;
+
+    libMesh::out << "\nProperty Values:\n" << properties_out.to_string() << std::endl;
+
+    libMesh::out << "Property Derivative w.r.t. DIM1\n" << dproperties_out.to_string() << std::endl;
+
+    return 0;
+}

src/base/CMakeLists.txt

@@ -45,7 +45,9 @@ target_sources(mast
         ${CMAKE_CURRENT_LIST_DIR}/transient_assembly.cpp
         ${CMAKE_CURRENT_LIST_DIR}/transient_assembly.h
         ${CMAKE_CURRENT_LIST_DIR}/transient_assembly_elem_operations.cpp
-        ${CMAKE_CURRENT_LIST_DIR}/transient_assembly_elem_operations.h)
+        ${CMAKE_CURRENT_LIST_DIR}/transient_assembly_elem_operations.h
+        ${CMAKE_CURRENT_LIST_DIR}/warping_assembly.cpp
+        ${CMAKE_CURRENT_LIST_DIR}/warping_assembly.h)
 
 configure_file(${CMAKE_CURRENT_LIST_DIR}/mast_config.h.in
                ${CMAKE_CURRENT_LIST_DIR}/mast_config.h)

src/base/boundary_condition_base.h

@@ -47,7 +47,8 @@ namespace MAST {
         EXHAUST,
         ISOTHERMAL,
         ADIABATIC,
-        BOUNDARY_VELOCITY
+        BOUNDARY_VELOCITY,
+        DOF_DIRICHLET
     };
 
 

src/base/field_function_base.h

@@ -120,6 +120,20 @@ namespace MAST {
             libmesh_error_msg("Must be implemented in derived class; " << __PRETTY_FUNCTION__ << " in " << __FILE__ << " at line " << __LINE__);
         }
 
+
+        /*!
+         *    calculates the value of the derivative of function with respect to
+         *    the function \p f at the specified point, \p p, and time,
+         *    \p t using finite differences, and returns it in \p v.
+         */
+        virtual void derivative (MAST::FunctionBase& f,
+                                 const libMesh::Point& p,
+                                 const Real t,
+                                 ValType& v) {
+
+            libmesh_error_msg("Must be implemented in derived class; " << __PRETTY_FUNCTION__ << " in " << __FILE__ << " at line " << __LINE__);
+        }
+
     protected:
 
     };

src/base/nonlinear_system.cpp

@@ -801,7 +801,9 @@ MAST::NonlinearSystem::adjoint_solve(const libMesh::NumericVector<Real>& X,
                                      MAST::AssemblyElemOperations&       elem_ops,
                                      MAST::OutputAssemblyElemOperations& output,
                                      MAST::AssemblyBase&                 assembly,
-                                     bool if_assemble_jacobian) {
+                                     bool if_assemble_jacobian,
+                                     bool compute_adjoint_rhs,
+                                     unsigned int i) {
 
 
     libmesh_assert(_operation == MAST::NonlinearSystem::NONE);
@@ -812,19 +814,21 @@ MAST::NonlinearSystem::adjoint_solve(const libMesh::NumericVector<Real>& X,
     LOG_SCOPE("adjoint_solve()", "NonlinearSystem");
 
     libMesh::NumericVector<Real>
-    &dsol  = this->add_adjoint_solution(),
-    &rhs   = this->add_adjoint_rhs();
+    &dsol  = this->add_adjoint_solution(i),
+    &rhs   = this->add_adjoint_rhs(i);
 
     assembly.set_elem_operation_object(elem_ops);
 
     if (if_assemble_jacobian)
         assembly.residual_and_jacobian(*solution, nullptr, matrix, *this);
 
-    assembly.calculate_output_derivative(X, if_localize_sol, output, rhs);
+    if (compute_adjoint_rhs)
+    {
+        assembly.calculate_output_derivative(X, if_localize_sol, output, rhs);
+        rhs.scale(-1.);
+    }
 
     assembly.clear_elem_operation_object();
-
-    rhs.scale(-1.);
 
     // Our iteration counts and residuals will be sums of the individual
     // results

src/base/nonlinear_system.h

@@ -157,7 +157,9 @@ namespace MAST {
                                    MAST::AssemblyElemOperations&       elem_ops,
                                    MAST::OutputAssemblyElemOperations& output,
                                    MAST::AssemblyBase&                 assembly,
-                                   bool if_assemble_jacobian           = true);
+                                   bool if_assemble_jacobian           = true,
+                                   bool compute_adjoint_rhs            = true,
+                                   unsigned int i                      = 0);
 
 
         /**

src/base/physics_discipline_base.cpp

@@ -23,6 +23,7 @@
 #include "base/parameter.h"
 #include "base/nonlinear_system.h"
 #include "boundary_condition/dirichlet_boundary_condition.h"
+#include "boundary_condition/dirichlet_dof_boundary_condition.h"
 #include "mesh/geom_elem.h"
 
 // libMesh includes
@@ -38,6 +39,7 @@ void
 MAST::PhysicsDisciplineBase::clear_loads() {
     _side_bc_map.clear();
     _vol_bc_map.clear();
+    _point_loads.clear();
 }
 
 
@@ -88,6 +90,13 @@ MAST::PhysicsDisciplineBase::add_dirichlet_bc(libMesh::boundary_id_type bid,
 }
 
 
+void
+MAST::PhysicsDisciplineBase::add_dirichlet_dof_bc(MAST::DOFDirichletBoundaryCondition& load)
+{
+    libmesh_assert(!_dirichlet_dof_bcs.count(&load));
+    _dirichlet_dof_bcs.insert(&load);
+}
+
 
 void
 MAST::PhysicsDisciplineBase::
@@ -229,6 +238,21 @@ init_system_dirichlet_bc(MAST::NonlinearSystem& sys) const {
 }
 
 
+void
+MAST::PhysicsDisciplineBase::
+init_system_dirichlet_dof_bc(MAST::NonlinearSystem& sys) {
+
+    // give the set of all dirichlet dof BCs to the MAST::DOFConstraint object 
+    // (a child class of libMesh::System:Constraint) 
+    _dof_constraint.reset(new MAST::DOFConstraint(_dirichlet_dof_bcs));
+
+    // tell the MAST::DOFConstraint which system it is acting on
+    _dof_constraint->setNonlinearSystem(sys);
+
+    // attach the MAST::DOFConstraint object to the system
+    sys.attach_constraint_object(*_dof_constraint);
+
+}
 
 
 void