Property Components
B2000++ property components provide the link between strain and stress or other element derivatives of an element. The property component will
Set up the proper property object for the material type.
Get the material data from the database MATERIAL table.
Provide the property->get_stress() function.
Elements will then request the element derivatives with the property->get_stress() method of the property component assigned to the element.
Input Processor
Unless special input definition is required, the input processor does not have to be modified. The MDL material block allows for specifying user-defined materials (type user “name”). Parameters are then (name type value) pairs. Example material:
material 1
type user USERMAT1 # case-sensitive!
key E F 71.e9 # key names must be uppercase!
key NU F 0.0
key DENSITY F 2700
key ALPHA F 21e-6
end
Property Implementation
Implementations of elements, element material properties, boundary
conditions, solvers, etc. are all sub-classes of b2000::Object.
Instances of b2000::ObjectType are factory objects which create
instances of specific sub-classes of b2000::Object and which, in
conjunction with b2000::Module, allow for the implementation of
plugins without having to recompile B2000++. b2000::Object and
b2000::ObjectType are explained in chapter Components.
The example implementation of the rod material property object is the
simplest possible: Unlike other B2000++ property implementations it gets
the material data straight from the RTable objects of the
MATERIAL table and from the ETAB.br.
Note
There may be some confusion about “material” and “property”.
“Element material” (from MATERIAL.mid dataset)
describes the material with constants, functions etc. “Property” in
the context of the B2000+++ code means a property object associated
to any element, providing the get_stress() function. The
property object can get information from MATERIAL.mid, PROPERTY.pid, and from
ETAB.br, depending on the element type.
The code starts with the include directives:
// The element property base class FortranElementProperty (inherits from
// ElementProperty). It makes the connection with factory and database and
// element. It
// (1) Creates *property* factory object,
// (2) loads data from the MATERIAL.<eid> dataset in the RTable object
// *rtable*,
// (3) gets the *element* attributes TL (cross-section area) and ECC
// (initial strain/stress/force) from ETAB.br.
#include "elements/properties/b2fortran_element_property.H"
#include "io/b2000_db/b2000_db_v3/b2domain_database.H"
//#include "utils/b2linear_algebra.H"
// The RTable object code
#include "utils/b2rtable.H"
// The base class of the property API for rod
// elements. ElementPropertyStressUSERMAT1 builds on SolidMechanicsRod
#include "elements/properties/b2solid_mechanics.H"
Declare the name space:
namespace b2000::b2dbv3 {
Now the actual code of the class
ElementPropertyStressUSERMAT1 starts:
//! User defined rod ElementProperty implementation example class.
class ElementPropertyStressUSERMAT1 :
virtual public FortranElementProperty,
public b2000::SolidMechanicsRod<double> // Is templatized.
// double impl.
{
The public section of ElementPropertyStressUSERMAT1. In
this implementation all public methods are embedded in the public
declaration:
public:
//! set_model_and_case() is called when the property object is
//! created. It loads the material data from the database.
void set_model_and_case(b2000::Model& model, b2000::Case& case_) {
#if VERBOSE > 0
std::cerr << "set_model_and_case of property "
<< material->get_object_name() << std::endl;
#endif
// rtable is the C++ representation of the user-defined
// material as specified in the MDL material block and stored
// on the database as a relational table MATERIAL.<eid>. This
// code block extracts the material data from the table and
// stores them in the private section of the object.
assert(material->number_of_layer() == 1); // Rod material
// cannot be layered
// cannot be layered
// ***Data from MATERIAL dataset***
const RTable& rtable = material->get_rtable(0); // layer 0
//const std::string sub_type = rtable.get_string("SUBTYPE");
// if (sub_type == "ISOTROPIC") {... // sub_type implicitly assumed
// to be ISOTROPIC!!!
E = rtable.get<double>("E");
ALPHA = rtable.get<double>("ALPHA", 0.);
DENSITY = rtable.get<double>("DENSITY", 0.);
C1 = rtable.get<double>("C1", 0.); // Constant of our
// "nonlinear" material
// **Data from ETAB (provided by the element code)***
// TL: cross section area for rod elements
if (tl.size() != 1)
Exception() << material->get_object_name() <<
": Properties size error, size found " << tl.size() <<
" requested: 1" << THROW;
if (tl[0] <= 0.0)
Exception() << material->get_object_name() <<
": Cross section area is <= 0" << THROW;
AREA = tl[0];
// Initial strain flag in ecc[0], initial strsin/stress/force
// in ecc[1]
if (ecc.size() != 0 && ecc[0] > 0)
INIT_STRAIN_FLAG = int(ecc[0]);
else
INIT_STRAIN_FLAG = 0;
#if VERBOSE > 0
std::cerr << " INIT_STRAIN_FLAG=" << INIT_STRAIN_FLAG << std::endl
<< " E=" << E << std::endl
<< " AREA=" << AREA << std::endl
<< " ALPHA=" << ALPHA<< std::endl
<< " DENSITY=" << DENSITY << std::endl
<< " C1=" << C1 << " (non-linear function test if != 0)"
<< std::endl;
#endif
} // end set_model_and_case
//! If the material model is linear or if it depends on the strain or
//! on the strain history.
SolidMechanics::LinearType linear(const int layer_id = -1) const {
return SolidMechanics::yes;
}
//! If the material is isotropic (check).
bool isotropic(const int layer_id = -1) const {
return true;
}
//! Given the current strain the function returns the stress
//! (force). In our example the force is returned. This is the
//! only method that the rod element calls.
void get_stress(
// in
b2000::Model* model,
const bool linear,
const EquilibriumSolution equilibrium_solution,
const double time,
const double delta_time,
b2000::GradientContainer* gradient_container,
b2000::SolverHints* solver_hints,
const Element* element,
const Vector<double, Vdense_constref> nodes_interpolation,
const double /*T*/ covariant_base[3],
const double /*T*/ displacement_gradient[3],
const double /*T*/ strain_xx,
const double /*T*/ velocity[3],
const double /*T*/ acceleration[3],
double /*T*/& force_x,
double /*T*/& C,
double /*T*/ inertia_force[3],
double /*T*/& density) {
#if VERBOSE > 0
std::cerr << " get_stress, element=" << element->get_id()
<< ", user material=" << material->get_object_name()
<< std::endl;
#endif
// Compute rod force force_x
C = E * AREA;
force_x = C * (strain_xx * (1.0 - C1 * strain_xx));
if (INIT_STRAIN_FLAG > 0) {
switch (int(ecc[0])) {
case 1:
force_x += C * ecc[1]; // initial strain
break;
case 2:
force_x += ecc[1] * AREA; // initial stress
break;
case 3:
force_x += ecc[1]; // inital force given
break;
default:
DBError() << "User material " << material->get_object_name() <<
", element " << element->get_id() <<
", Invalid initial stress/strain/force type " <<
int(ecc[0].real()) << THROW;
break;
}
}
#if VERBOSE > 1
std::cerr << " C=E*AREA=" << C << ", force_x=" << force_x << std::endl;
#endif
double temp = 0;
{
double d_temperature_d_time = 0;
const double one = 1;
get_material_temperature(
element, nodes_interpolation,
1, &one, time, temp, d_temperature_d_time);
}
force_x -= C * ALPHA * temp;
if (inertia_force) {
std::fill_n(inertia_force, 3, 0.0);
if (acceleration) {
double densitya = DENSITY * AREA;
inertia_force[0] = densitya * acceleration[0];
inertia_force[1] = densitya * acceleration[1];
inertia_force[2] = densitya * acceleration[2];
}
#if VERBOSE > 1
std::cerr << " acceleration=(" << acceleration[0] << ","
<< acceleration[1] << ","
<< acceleration[2] << ")" << std::endl;
std::cerr << " inertia forces=(" << inertia_force[0] << ","
<< inertia_force[1] << ","
<< inertia_force[2] << ")" << std::endl;
#endif
}
// If no gradients requested or during iteration
if (gradient_container == 0)
return;
#if VERBOSE > 1
std::cerr << " store gradients " << std::endl;
#endif
//
// Storage to gradient container
//
double /*T*/ s; // cauchy_stress conversion
if (linear)
s = 1;
else {
double /*T*/ tmp1 = 0;
double /*T*/ tmp2 = 0;
for (int i = 0; i != 3; ++i) {
tmp1 += covariant_base[i] * covariant_base[i];
const double /*T*/ tmp3 = covariant_base[i] + displacement_gradient[i];
tmp2 += tmp3 * tmp3;
}
s = std::sqrt(tmp2 / tmp1);
}
// Store the strain eps_xx
const GradientContainer::FieldDescription strain_descr = {
"STRAIN_SECTION_ROD",
"Exx",
(linear ? "Linear strain in the element reference frame" :
"Green Lagrange strain in the element reference frame"),
3, 1, -1, -1, false, typeid(double /*T*/)
};
if (gradient_container->compute_field_value(strain_descr.name))
gradient_container->set_field_value(strain_descr, &strain_xx);
// Store the stress sigma_xx
const GradientContainer::FieldDescription stress_descr = {
"STRESS_SECTION_ROD",
"Sxx",
(linear ? "Linear stress in the element reference frame" :
"Cauchy stress in the element reference frame"),
3, 1, -1, -1, false, typeid(double /*T*/)
};
if (gradient_container->compute_field_value(stress_descr.name)) {
double /*T*/ cauchy_stress_xx = s * force_x / AREA;
gradient_container->set_field_value(stress_descr, &cauchy_stress_xx);
}
// Store the rod force
const GradientContainer::FieldDescription force_descr = {
"FORCE_ELEMENT",
"Fx",
("Section force"),
3, 1, -1, -1, false, typeid(double /*T*/)
};
if (gradient_container->compute_field_value(force_descr.name)) {
double /*T*/ cauchy_force_x = s * force_x;
gradient_container->set_field_value(force_descr, &cauchy_force_x);
}
density = get_density(); // return density
} // end get_stress
double /*T*/ get_area() const {
return AREA; // UnimplementedError() << THROW;
//return 0;
}
double /*T*/ get_density() const {
return DENSITY;
}
typedef ObjectTypeComplete<
ElementPropertyStressUSERMAT1,
FortranElementProperty::type_t> type_t;
static type_t type;
The private section of ElementPropertyStressUSERMAT1. In this
implementation only material data are stored (no methods).
private:
// These are extracted from Rtable in set_model_and_case()
double AREA;
double E;
double ALPHA;
double DENSITY;
double C1;
int INIT_STRAIN_FLAG;
}; // end ElementPropertyStressUSERMAT1
// Registration of the class under its name. Thus, given the material
// type name, b2000++ can instantiate the element objects. This must
// match the name of the shared library (see CMakeLists.txt).
Finally the object declared globally for the element factory:
ElementPropertyStressUSERMAT1::type_t
ElementPropertyStressUSERMAT1::type(
"ElementPropertyStressUSERMAT1",
"", StringList("USERMAT1"),
element::module, &b2000::ElementProperty::type);
} // namespace b2000::b2dbv3
Build the Property Project
The user property is built with the cmake build system. The output of the build process is the shared library
<CMAKE_INSTALL_PREFIX>/lib/libb2000.element.USERMAT1.so
To build the project, go to the components/rod_property
directory, and
cd <b2programming>/components/rod_property
mkdir -p _build
cd _build
cnake [-DCMAKE_INSTALL_PREFIX=$SMR_PREFIX] ..
make install
The example cmake user written property project found in the
b2programming package under directory
components/rod_property is set up as follows:
rod_property: CMakeLists.txt
rod_property/src: b2element_property_stress_rod_user.C, CMakeLists.txt
rod_property/test: rod45, nonlinear_elastic
Root cmake file rod_property/CMakeList.txt:
# cmake file to build user-defined property
# =========================================
# Example configuration under this directory. Note: The target "make
# install" must be executed.
# Default cmake with $SMR_PREFIX and $SMR_CMAKEMODULE defined:
# mkdir _build; cd _build; cmake ..; make install
# Options to cmake:
# -DCMAKE_INSTALL_PREFIX=<path-to-install-prefix>
# -DCMAKE_MODULE_PATH=<path-to-cmake-modules>
# Author: S. R. P. Merazzi, SMR SA
# Project
# -------
cmake_minimum_required(VERSION 3.22)
project(property DESCRIPTION "user rod property" VERSION 1.0.0 LANGUAGES CXX)
include("../common.cmake")
add_subdirectory(src)
cmake file rod_element/src/CMakeList.txt for
building the library:
# cmake file to generate B2000++ make file for user-written rod property
# ======================================================================
# Authors: Silvio Merazzi, SMR SA, Thomas Rothermel, DLR
# Create a shared library for the new USER1 material. It will be
# loaded dynamically by B2000++ solver if a USER1 property is
# found. Compile and add b2element_cable_user.C.
# Create the new library
add_library(b2000.element.USERMAT1 SHARED "")
# Compile and add to library
target_sources(b2000.element.USERMAT1
PRIVATE
b2element_property_stress_rod_user.C
)
# Link library with B2000 and Python libraries
target_link_libraries(b2000.element.USERMAT1
PUBLIC B2000::B2000
PUBLIC Python3::Python
)
# Installs the user element library
Install(
TARGETS b2000.element.USERMAT1
DESTINATION ${CMAKE_INSTALL_LIBDIR}
)
And the project $B2PROGRAMMING/components/rod_element/CMakeList.txt file:
#[[
Main cmake file to build user-defined element
=============================================
Example configuration with default $SMR_PREFIX and $SMR_CMAKEMODULE defined:
mkdir _build; cd _build; cmake ..; make install
Note: The target "make install" must be executed.
Options to cmake:
-DCMAKE_INSTALL_PREFIX=<path-to-install-prefix>
-DCMAKE_MODULE_PATH=<path-to-cmodules>
Author: S. R. P. Merazzi, SMR SA
]]
cmake_minimum_required(VERSION 3.19)
project(rod DESCRIPTION "user-defined rod element" VERSION 1.0.0 LANGUAGES CXX)
include("../common.cmake")
add_subdirectory(src)
Testing the Property Implementation
If the build is successful, select one of the tests found under
rod_property/test and run the test with your b2000++
application.
Test case “rod45”
Example with two rods loaded in compression by a force. The problem is
a pure 2D problem, with all DOFs except Ux and Uy removed. Note that
the test is run with the user element R2.S.USER: The rod_element
project has to be built prior to executing this test.
The user defined material parameters are specified in the
rod.mdl:
material 1
type user USERMAT1 # case-sensitive!
key E F 71.e9 # key names must be uppercase!
key NU F 0.0
key DENSITY F 2700
key ALPHA F 21e-6
end
They are identical to the built-in isotropic material specification of
material 2, see rod.mdl:
material 2
type isotropic
E 71.e9
NU 0.0
DENSITY 2700
ALPHA 21e-6
end
make
Execute with user material 1:
b2000++ rod.mdl
b2print_solutions rod.b2m
or
make
Should produce the output
DOF field DISP, branch=1, case=1 cycle=0, nodeset="all nodes"
EID Ux Uy Uz
1 0 0 0
2 0 -0.01408 0
3 0 0 0
Expected Uy(node 2)=-0.01408
Setting the element material identifier “mid” to 2 in rod.mdl
and running the test will give the same results. Try out.
Test case “rod10”)
The non-linear system test rod_property/test/rod10
demonstrates the ‘snap-through’ effect (see plot below). The model is
similar to the rod45 model, with a smaller ration of height to with,
provoking the ‘snap-through’.
rod10: Model.
The USERMAT1 user material has an assumed non-linear (hypoelastic)
behavior
built in the function get_stress of the USERMAT1 material. C++
code ‘snippet’:
// Compute rod force force_x
C = E * AREA;
force_x = C * (strain_xx * (1.0 - C1 * strain_xx));
Note that we work with the rod axial force instead of the stress. And if \(C_1\) is set to 0 the material is linear elastic.
With the tension test FE model tension.mdl the stress in the rods
must be equal to the stress-strain (”\(\sigma_{xx} = E
(\epsilon_{xx}\)”) curve
Element stress-strain curve (made by the Simples script
make-strain-stress-graph.py).