diff --git a/.cmake/mito_tests_mito_lib.cmake b/.cmake/mito_tests_mito_lib.cmake index 88ab480e6..b771cb58d 100644 --- a/.cmake/mito_tests_mito_lib.cmake +++ b/.cmake/mito_tests_mito_lib.cmake @@ -46,11 +46,15 @@ mito_test_driver(tests/mito.lib/discrete/mesh_field.cc) # fem mito_test_driver(tests/mito.lib/fem/block_grad_grad.cc) +mito_test_driver(tests/mito.lib/fem/block_grad_grad_segment.cc) mito_test_driver(tests/mito.lib/fem/block_mass.cc) +mito_test_driver(tests/mito.lib/fem/block_mass_segment.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_construction.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p1.cc) mito_test_driver(tests/mito.lib/fem/shape_functions_triangle_p2.cc) mito_test_driver(tests/mito.lib/fem/isoparametric_triangle.cc) +mito_test_driver(tests/mito.lib/fem/shape_functions_segment_p1.cc) +mito_test_driver(tests/mito.lib/fem/isoparametric_segment.cc) # io mito_test_driver(tests/mito.lib/io/summit_mesh_reader_2D.cc) diff --git a/AUTHORS b/AUTHORS index e1694c831..76a562312 100644 --- a/AUTHORS +++ b/AUTHORS @@ -2,3 +2,4 @@ This is a list of the MiTo copyright holders: Bianca Giovanardi Anwar Koshakji Michael A.G. Aïvázis + Pawel J. Rok diff --git a/lib/mito/fem/elements/IsoparametricSegment.h b/lib/mito/fem/elements/IsoparametricSegment.h new file mode 100644 index 000000000..3c867e35d --- /dev/null +++ b/lib/mito/fem/elements/IsoparametricSegment.h @@ -0,0 +1,74 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// DESIGN NOTES +// Class {IsoparametricSegment} represents a first order simplex (segment) equipped with parametric coordinates. + + +namespace mito::fem { + + class IsoparametricSegment : public utilities::Invalidatable { + public: + // the dimension of the physical space + static constexpr int dim = 1; + // the discretization node type + using discretization_node_type = discrete::discretization_node_t; + // the underlying cell type + using cell_type = geometry::segment_t; + + protected: + // cartesian coordinates in 1D + using coordinates_type = geometry::coordinates_t<1, geometry::CARTESIAN>; + // the coordinate system type + using coordinate_system_type = geometry::coordinate_system_t; + // the vector type + using vector_type = tensor::vector_t<1>; + + public: + // the default constructor + constexpr IsoparametricSegment( + const cell_type & cell, const coordinate_system_type & coord_system) : + _cell(cell), + _x0{ coord_system.coordinates(cell.nodes()[0]->point()) - coordinates_type{} }, + _x1{ coord_system.coordinates(cell.nodes()[1]->point()) - coordinates_type{} } + {} + + // destructor + constexpr ~IsoparametricSegment() = default; + + // delete move constructor + constexpr IsoparametricSegment(IsoparametricSegment &&) noexcept = delete; + + // delete copy constructor + constexpr IsoparametricSegment(const IsoparametricSegment &) = delete; + + // delete assignment operator + constexpr IsoparametricSegment & operator=(const IsoparametricSegment &) = delete; + + // delete move assignment operator + constexpr IsoparametricSegment & operator=(IsoparametricSegment &&) noexcept = delete; + + public: + // get the geometric simplex + constexpr auto cell() const noexcept -> const cell_type & { return _cell; } + + protected: + // QUESTION: do we need to maintain a reference to the geometric simplex? + // a const reference to the geometric simplex + const cell_type & _cell; + + // the coordinates of the discretization nodes of the segment + const vector_type _x0; + const vector_type _x1; + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/elements/elements_library.h b/lib/mito/fem/elements/elements_library.h index d52cfc68d..ac2cb50f4 100644 --- a/lib/mito/fem/elements/elements_library.h +++ b/lib/mito/fem/elements/elements_library.h @@ -7,6 +7,7 @@ #pragma once +#include "seg1/public.h" #include "tri1/public.h" #include "tri2/public.h" diff --git a/lib/mito/fem/elements/public.h b/lib/mito/fem/elements/public.h index 40f144579..1b33bd6be 100644 --- a/lib/mito/fem/elements/public.h +++ b/lib/mito/fem/elements/public.h @@ -17,6 +17,7 @@ #include "api.h" // classes implementation +#include "IsoparametricSegment.h" #include "IsoparametricTriangle.h" #include "Discretizer.h" diff --git a/lib/mito/fem/elements/seg1/DiscretizerCG.h b/lib/mito/fem/elements/seg1/DiscretizerCG.h new file mode 100644 index 000000000..0735af4c1 --- /dev/null +++ b/lib/mito/fem/elements/seg1/DiscretizerCG.h @@ -0,0 +1,65 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::fem { + + // discretizer specialization for {IsoparametricSegmentP1} with continuous Galerkin + template <> + struct Discretizer { + template < + typename manifoldT, typename constraintsT, typename elements_type, typename map_type, + typename constrained_nodes_type> + static void apply( + const manifoldT & manifold, const constraintsT & constraints, elements_type & elements, + map_type & node_map, constrained_nodes_type & constrained_nodes) + { + + // the discretization node type + using discretization_node_type = + typename IsoparametricSegmentP1::discretization_node_type; + + // the connectivity type + using connectivity_type = typename IsoparametricSegmentP1::connectivity_type; + + // get the coordinate system of the manifold + const auto & coord_system = manifold.coordinate_system(); + + // loop on the cells of the mesh + for (const auto & cell : manifold.elements()) { + + // get the nodes of the cell + const auto & nodes = cell.nodes(); + + // add the nodes to the map (if the mesh node is already present in the map, + // then the present discretization node is used) + auto node_0 = + node_map.insert({ nodes[0], discretization_node_type() }).first->second; + auto node_1 = + node_map.insert({ nodes[1], discretization_node_type() }).first->second; + + // create a finite element for each cell and add it to the pile + elements.emplace(cell, coord_system, connectivity_type{ node_0, node_1 }); + } + + // populate the constrained nodes + // In 1D, constraints.domain() is a set of nodes, not a mesh with cells, so we loop on the nodes directly + for (const auto & node : constraints.domain()) { + // get the discretization node associated with the mesh node from the map + auto it = node_map.find(node); + // add the node to the constrained nodes + constrained_nodes.insert(it->second); + } + + // all done + return; + } + }; +} + +// end of file diff --git a/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h new file mode 100644 index 000000000..c4e8dde6b --- /dev/null +++ b/lib/mito/fem/elements/seg1/IsoparametricSegmentP1.h @@ -0,0 +1,136 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// DESIGN NOTES +// Class {IsoparametricSegmentP1} represents a first order simplex (segment) living in 1D cartesian space, +// equipped with linear shape functions defined in the parametric space. + + +namespace mito::fem { + + class IsoparametricSegmentP1 : public IsoparametricSegment { + + public: + // the degree of the finite element + static constexpr int degree = 1; + // the type of shape functions + using shape_functions_type = ShapeSegmentP1; + // the canonical element type + using canonical_element_type = typename shape_functions_type::reference_element_type; + // the parametric coordinates type + using parametric_coordinates_type = + typename canonical_element_type::parametric_coordinates_type; + // the linear shape functions + static constexpr auto shape_functions = shape_functions_type(); + // the number of discretization nodes + static constexpr int n_nodes = shape_functions_type::N; + // a collection of discretization nodes + using connectivity_type = std::array; + + public: + // the default constructor + inline IsoparametricSegmentP1( + const cell_type & geometric_simplex, const coordinate_system_type & coord_system, + const connectivity_type & connectivity) : + IsoparametricSegment(geometric_simplex, coord_system), + _connectivity(connectivity) + {} + + // destructor + inline ~IsoparametricSegmentP1() = default; + + // delete move constructor + constexpr IsoparametricSegmentP1(IsoparametricSegmentP1 &&) noexcept = delete; + + // delete copy constructor + constexpr IsoparametricSegmentP1(const IsoparametricSegmentP1 &) = delete; + + // delete assignment operator + constexpr IsoparametricSegmentP1 & operator=(const IsoparametricSegmentP1 &) = delete; + + // delete move assignment operator + constexpr IsoparametricSegmentP1 & operator=(IsoparametricSegmentP1 &&) noexcept = delete; + + public: + // get the discretization nodes + constexpr auto connectivity() const noexcept -> const connectivity_type & + { + return _connectivity; + } + + // get the isoparametric mapping from parametric coordinates to physical coordinates + constexpr auto parametrization() const + { + // get the shape functions + constexpr auto phi_0 = shape_functions.shape<0>(); + constexpr auto phi_1 = shape_functions.shape<1>(); + + // return the isoparametric mapping from parametric to physical coordinates + return mito::fields::linear_combination( + std::array{ _x0, _x1 }, phi_0, phi_1); + } + + // get the shape function associated with local node {a} + template + requires(a >= 0 && a < n_nodes) + constexpr auto shape() const + { + // return the shape functions + return shape_functions.shape(); + } + + // get the jacobian of the isoparametric mapping from parametric to actual coordinates + constexpr auto jacobian() const + { + // assemble the jacobian as a function of parametric coordinates + auto jacobian_function = functions::function( + [&](const parametric_coordinates_type & xi) -> tensor::matrix_t<1> { + // get the shape functions derivatives + constexpr auto dphi_0 = shape_functions.dshape<0>(); + constexpr auto dphi_1 = shape_functions.dshape<1>(); + + // compute the jacobian of the isoparametric mapping: dx/dxi + auto dx_dxi = _x0 * dphi_0(xi) + _x1 * dphi_1(xi); + // wrap the result in a 1x1 matrix + return tensor::matrix_t<1>{ dx_dxi }; + }); + + // and return it + return jacobian_function; + } + + // get the gradient of the a-th shape function as a function of parametric coordinates + template + requires(a >= 0 && a < n_nodes) + constexpr auto gradient() const + { + // assemble the gradient as a function of parametric coordinates + auto gradient_function = + fields::field([&](const parametric_coordinates_type & xi) -> tensor::vector_t<1> { + // the jacobian of the mapping from the reference element to the physical + // element evaluated at {xi} + auto J = jacobian()(xi); + // the derivative of the coordinates with respect to the parametric coordinates + auto J_inv = tensor::inverse(J); + // return the spatial gradients of the shape functions evaluated at {xi} + return shape_functions.dshape()(xi) * J_inv; + }); + // and return it + return gradient_function; + } + + private: + // the discretization nodes of the simplex + const connectivity_type _connectivity; + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/elements/seg1/ShapeSegmentP1.h b/lib/mito/fem/elements/seg1/ShapeSegmentP1.h new file mode 100644 index 000000000..c4c576a54 --- /dev/null +++ b/lib/mito/fem/elements/seg1/ShapeSegmentP1.h @@ -0,0 +1,58 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::fem { + + class ShapeSegmentP1 { + + public: + // the reference element + using reference_element_type = geometry::reference_segment_t; + // the number of shape functions + static constexpr int N = 2; + + private: + // linear shape functions on the reference segment in parametric coordinates + static constexpr auto xi_0 = reference_element_type::xi<0>; + + // linear shape functions on the segment + static constexpr auto phi_0 = 1.0 - xi_0; + static constexpr auto phi_1 = xi_0; + + // the shape functions + static constexpr auto phi = std::make_tuple(phi_0, phi_1); + + // the gradients of the shape functions + static constexpr auto dphi = std::make_tuple( + fields::gradient(phi_0), fields::gradient(phi_1)); + + public: + // get the a-th shape function as a function of parametric coordinates + template + requires(a >= 0 && a < N) + constexpr auto shape() const + { + // return the a-th shape function + return std::get(phi); + } + + // get the a-th shape function's derivative as a function of parametric coordinates + template + requires(a >= 0 && a < N) + constexpr auto dshape() const + { + // return the a-th shape function's gradient + return std::get(dphi); + } + }; + +} // namespace mito + + +// end of file diff --git a/lib/mito/fem/elements/seg1/api.h b/lib/mito/fem/elements/seg1/api.h new file mode 100644 index 000000000..3fcb2b63b --- /dev/null +++ b/lib/mito/fem/elements/seg1/api.h @@ -0,0 +1,21 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +namespace mito::fem { + + // specialization for linear shape functions on segments in 1D + template <> + struct isoparametric_simplex<1, geometry::segment_t<1>> { + using type = IsoparametricSegmentP1; + }; + +} + + +// end of file diff --git a/lib/mito/fem/elements/seg1/public.h b/lib/mito/fem/elements/seg1/public.h new file mode 100644 index 000000000..d3f42ddbf --- /dev/null +++ b/lib/mito/fem/elements/seg1/public.h @@ -0,0 +1,19 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +// code guard +#pragma once + + +// classes implementation +#include "ShapeSegmentP1.h" +#include "IsoparametricSegmentP1.h" +#include "DiscretizerCG.h" + +// published types and factories +#include "api.h" + + +// end of file diff --git a/tests/mito.lib/fem/block_grad_grad_segment.cc b/tests/mito.lib/fem/block_grad_grad_segment.cc new file mode 100644 index 000000000..9e25f1c2b --- /dev/null +++ b/tests/mito.lib/fem/block_grad_grad_segment.cc @@ -0,0 +1,72 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates +using coordinates_t = mito::geometry::coordinates_t<1, mito::geometry::CARTESIAN>; +// the type of coordinate system +using coord_system_t = mito::geometry::coordinate_system_t; +// the type of discretization node +using discretization_node_t = mito::discrete::discretization_node_t; +// the type of cell +using cell_t = mito::geometry::segment_t<1>; +// the reference simplex +using reference_simplex_t = mito::geometry::reference_segment_t; +// Gauss quadrature on segments with degree of exactness 2 +using quadrature_rule_t = + mito::quadrature::quadrature_rule_t; + +// instantiate the quadrature rule +constexpr auto quadrature_rule = quadrature_rule_t(); + + +TEST(Fem, BlockGradGradSegment) +{ + // the coordinate system + auto coord_system = coord_system_t(); + + // build nodes + auto node_0 = mito::geometry::node(coord_system, { 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<1>({ node_0, node_1 }); + + { + // first order isoparametric segment + using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>; + + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + + // a finite element + auto element_p1 = element_p1_t( + geometric_simplex, coord_system, + { discretization_node_0, discretization_node_1 }); + + // a grad-grad matrix block + auto grad_grad_block = + mito::fem::blocks::grad_grad_block(); + + // the analytical elementary stiffness matrix + auto analytical_block = mito::tensor::matrix_t<2>{ 1.0, -1.0, -1.0, 1.0 }; + + // compute the elementary contribution of the block + auto computed_block = grad_grad_block.compute(element_p1); + + // compute the error + auto error = mito::tensor::norm(computed_block - analytical_block); + + // check the error is zero to machine precision + EXPECT_DOUBLE_EQ(0.0, error); + } + + // all done + return; +} \ No newline at end of file diff --git a/tests/mito.lib/fem/block_mass_segment.cc b/tests/mito.lib/fem/block_mass_segment.cc new file mode 100644 index 000000000..97224cb00 --- /dev/null +++ b/tests/mito.lib/fem/block_mass_segment.cc @@ -0,0 +1,72 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates +using coordinates_t = mito::geometry::coordinates_t<1, mito::geometry::CARTESIAN>; +// the type of coordinate system +using coord_system_t = mito::geometry::coordinate_system_t; +// the type of discretization node +using discretization_node_t = mito::discrete::discretization_node_t; +// the type of cell +using cell_t = mito::geometry::segment_t<1>; +// the reference simplex +using reference_simplex_t = mito::geometry::reference_segment_t; +// Gauss quadrature on segments with degree of exactness 2 +using quadrature_rule_t = + mito::quadrature::quadrature_rule_t; + +// instantiate the quadrature rule +constexpr auto quadrature_rule = quadrature_rule_t(); + + +TEST(Fem, BlockMassSegment) +{ + // the coordinate system + auto coord_system = coord_system_t(); + + // build nodes + auto node_0 = mito::geometry::node(coord_system, { 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<1>({ node_0, node_1 }); + + { + // first order isoparametric segment + using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>; + + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + + // a finite element + auto element_p1 = element_p1_t( + geometric_simplex, coord_system, + { discretization_node_0, discretization_node_1 }); + + // a mass matrix block + auto mass_block = mito::fem::blocks::mass_block(); + + // the analytical elementary mass matrix + auto analytical_block = + 1.0 / 6.0 * mito::tensor::matrix_t<2>{ 2.0, 1.0, 1.0, 2.0 }; + + // compute the elementary contribution of the block + auto computed_block = mass_block.compute(element_p1); + + // compute the error + auto error = mito::tensor::norm(computed_block - analytical_block); + + // check the error is reasonable + EXPECT_NEAR(0.0, error, 1.5e-16); + } + + // all done + return; +} diff --git a/tests/mito.lib/fem/isoparametric_segment.cc b/tests/mito.lib/fem/isoparametric_segment.cc new file mode 100644 index 000000000..d5f5a6ba2 --- /dev/null +++ b/tests/mito.lib/fem/isoparametric_segment.cc @@ -0,0 +1,124 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// the type of coordinates +using coordinates_t = mito::geometry::coordinates_t<1, mito::geometry::CARTESIAN>; +// the type of coordinate system +using coord_system_t = mito::geometry::coordinate_system_t; +// the type of discretization node +using discretization_node_t = mito::discrete::discretization_node_t; +// the type of cell +using cell_t = mito::geometry::segment_t<1>; +// the reference simplex +using reference_simplex_t = mito::geometry::reference_segment_t; +// Gauss quadrature on segments with degree of exactness 2 +using quadrature_rule_t = + mito::quadrature::quadrature_rule_t; + + +// instantiate the quadrature rule +constexpr auto quadrature_rule = quadrature_rule_t(); + + +// test that all shape functions sum to 1.0 at any quadrature point +auto +test_partition_of_unity(const auto & element) +{ + // the number of quadrature points per element + constexpr int n_quads = quadrature_rule_t::npoints; + + // the number of nodes per element + constexpr int n_nodes = mito::utilities::base_type::n_nodes; + + // loop on the quadrature points + mito::tensor::constexpr_for_1([&]() { + // the parametric coordinates of the quadrature point + constexpr auto xi = quadrature_rule.point(q); + + // compute the sum of the shape functions at {xi} for all nodes + constexpr auto sum = + ([]( + const auto & element, const auto & xi, mito::tensor::integer_sequence) { + return ((element.template shape()(xi)) + ...); + })(element, xi, mito::tensor::make_integer_sequence{}); + + // check the sum of the shape functions + static_assert(1.0 == sum); + }); + + // all done + return; +} + +// test that the gradients of all shape functions sum to 0.0 at any quadrature point +auto +test_gradient_consistency(const auto & element) +{ + // the number of quadrature points per element + constexpr int n_quads = quadrature_rule_t::npoints; + + // the number of nodes per element + constexpr int n_nodes = mito::utilities::base_type::n_nodes; + + // loop on the quadrature points + mito::tensor::constexpr_for_1([&]() { + // the parametric coordinates of the quadrature point + constexpr auto xi = quadrature_rule.point(q); + + // compute the sum of the shape functions gradients at {xi} for all nodes + auto sum = + ([]( + const auto & element, const auto & xi, mito::tensor::integer_sequence) { + return ((element.template gradient()(xi)) + ...); + })(element, xi, mito::tensor::make_integer_sequence{}); + + // check the sum of the shape functions gradients + EXPECT_NEAR(0.0, sum, 3.0e-16); + }); + + // all done + return; +} + + +TEST(Fem, IsoparametricSegment) +{ + // the coordinate system + auto coord_system = coord_system_t(); + + // build nodes + auto node_0 = mito::geometry::node(coord_system, { 0.0 }); + auto node_1 = mito::geometry::node(coord_system, { 1.0 }); + + // make a geometric simplex + auto geometric_simplex = mito::geometry::segment<1>({ node_0, node_1 }); + + { + // first order isoparametric segment + using element_p1_t = mito::fem::isoparametric_simplex_t<1, cell_t>; + + // build the discretization nodes + auto discretization_node_0 = discretization_node_t(); + auto discretization_node_1 = discretization_node_t(); + + // a finite element + auto element_p1 = element_p1_t( + geometric_simplex, coord_system, + { discretization_node_0, discretization_node_1 }); + + // check that first order shape functions are a partition of unity + test_partition_of_unity(element_p1); + + // check that the gradients of first order shape functions sum to 0.0 + test_gradient_consistency(element_p1); + } + + // all done + return; +} diff --git a/tests/mito.lib/fem/shape_functions_segment_p1.cc b/tests/mito.lib/fem/shape_functions_segment_p1.cc new file mode 100644 index 000000000..f3388574a --- /dev/null +++ b/tests/mito.lib/fem/shape_functions_segment_p1.cc @@ -0,0 +1,42 @@ +// -*- c++ -*- +// +// Copyright (c) 2020-2024, the MiTo Authors, all rights reserved +// + +#include +#include + + +// first order shape functions type +using shape_t = mito::fem::ShapeSegmentP1; +// the parametric coordinates type +using parametric_coordinates_type = shape_t::reference_element_type::parametric_coordinates_type; + + +TEST(Fem, ShapeSegmentP1) +{ + // first order shape functions + constexpr auto element = shape_t(); + + // node 0 in parametric coordinates + constexpr auto n0 = parametric_coordinates_type{ 0.0 }; + // node 1 in parametric coordinates + constexpr auto n1 = parametric_coordinates_type{ 1.0 }; + + // the shape function associated with local node {0} + constexpr auto phi_0 = element.shape<0>(); + // check that the shape function at node 0 is 1.0 + static_assert(1.0 == phi_0(n0)); + // check that the shape function at node 1 is 0.0 + static_assert(0.0 == phi_0(n1)); + + // the shape functions at node 1 + constexpr auto phi_1 = element.shape<1>(); + // check that the shape function at node 0 is 0.0 + static_assert(0.0 == phi_1(n0)); + // check that the shape function at node 1 is 1.0 + static_assert(1.0 == phi_1(n1)); + + // all done + return; +} diff --git a/tests/mito.lib/fem/verification_segment_p1.py b/tests/mito.lib/fem/verification_segment_p1.py new file mode 100644 index 000000000..de134cd19 --- /dev/null +++ b/tests/mito.lib/fem/verification_segment_p1.py @@ -0,0 +1,38 @@ +import sympy +from sympy.vector import CoordSys3D, gradient + +# define the coordinate system +C = CoordSys3D('C') + +# define the coordinate +x = C.x + +# define the shape functions +phi0 = 1 - x +phi1 = x + +# report +print("Computing the mass matrix for the 2 shape functions of a segment element:") + +# compute the mass matrix via double integration +phis = [phi0, phi1] +print("Mass matrix:") +for i, phi_i in enumerate(phis): + for j, phi_j in enumerate(phis): + integral = sympy.integrate(phi_i * phi_j, (x, 0, 1)) + print(f"M[{i},{j}] = {integral}") + +# report +print("Computing the stiffness matrix for the 2 shape functions of a segment element:") + +# compute the shape functions gradients +dphi0 = gradient(phi0) +dphi1 = gradient(phi1) + +# compute the stiffness matrix via double integration +dphis = [dphi0, dphi1] +print("Stiffness matrix:") +for i, dphi_i in enumerate(dphis): + for j, dphi_j in enumerate(dphis): + integral = sympy.integrate(dphi_i & dphi_j, (x, 0, 1)) + print(f"K[{i},{j}] = {integral}")