diff --git a/CMakeLists.txt b/CMakeLists.txt index 392354d..01c0c21 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -17,9 +17,9 @@ endif() include(gsConfig) ## Collect files -aux_header_directory (${CMAKE_CURRENT_SOURCE_DIR} ${PROJECT_NAME}_H ) -aux_cpp_directory (${CMAKE_CURRENT_SOURCE_DIR} ${PROJECT_NAME}_CPP) -aux_tmpl_header_directory(${CMAKE_CURRENT_SOURCE_DIR} ${PROJECT_NAME}_HPP) +aux_header_directory (${CMAKE_CURRENT_SOURCE_DIR}/src ${PROJECT_NAME}_H ) +aux_cpp_directory (${CMAKE_CURRENT_SOURCE_DIR}/src ${PROJECT_NAME}_CPP) +aux_tmpl_header_directory(${CMAKE_CURRENT_SOURCE_DIR}/src ${PROJECT_NAME}_HPP) if( (NOT GISMO_BUILD_LIB) ) aux_instance_directory (${CMAKE_CURRENT_SOURCE_DIR} ${PROJECT_NAME}_INS) @@ -45,6 +45,10 @@ set_target_properties(${PROJECT_NAME} PROPERTIES set(gismo_MODULES ${gismo_MODULES} $ CACHE INTERNAL "G+Smo modules" ) +#Symlink include dir (in case your headers are in /src) +execute_process( COMMAND ${CMAKE_COMMAND} -E create_symlink ${CMAKE_CURRENT_SOURCE_DIR}/src ${CMAKE_BINARY_DIR}/${PROJECT_NAME}) + + install(DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}/${PROJECT_NAME}" DESTINATION include/gismo FILES_MATCHING PATTERN "*.h" ) diff --git a/examples/beam_nonLinElast_3D.cpp b/examples/beam_nonLinElast_3D.cpp new file mode 100644 index 0000000..0233841 --- /dev/null +++ b/examples/beam_nonLinElast_3D.cpp @@ -0,0 +1,125 @@ +/// This is an example of using the nonlinear elasticity solver on a 3D multi-patch geometry +/// The problems is part of the EU project "Terrific". +/// +/// Authors: O. Weeger (2012-1015, TU Kaiserslautern), +/// A.Shamanskiy (2016 - ...., TU Kaiserslautern) +#include +#include +#include +#include +#include +#include + +using namespace gismo; + +int main(int argc, char* argv[]){ + + gsInfo << "Testing the nonlinear elasticity solver in 3D.\n"; + + //=====================================// + // Input // + //=====================================// + + real_t youngsModulus = 74e9; + real_t poissonsRatio = 0.33; + index_t materialLaw = material_law::saint_venant_kirchhoff; + index_t numUniRef = 0; + index_t numDegElev = 0; + index_t numPlotPoints = 10000; + + // minimalistic user interface for terminal + gsCmdLine cmd("Testing the linear elasticity solver in 3D."); + cmd.addInt("l","law","Material law: 0 - St.V.-K., 1 - neoHookeLn, 2 - neoHookeQuad",materialLaw); + cmd.addInt("r","refine","Number of uniform refinement application",numUniRef); + cmd.addInt("d","degelev","Number of degree elevation application",numDegElev); + cmd.addInt("p","points","Number of points to plot to Paraview",numPlotPoints); + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + //=============================================// + // Scanning geometry and creating bases // + //=============================================// + + // scanning geometry + gsMultiPatch<> patch, geometry; + patch.addPatch(gsNurbsCreator<>::BSplineCube()); + patch.patch(0).coefs().col(0)*=10; + geometry = patch.uniformSplit(0); + geometry.computeTopology(); + // creating basis + gsMultiBasis<> basis(geometry); + for (index_t i = 0; i < numDegElev; ++i) + basis.degreeElevate(); + for (index_t i = 0; i < numUniRef; ++i) + basis.uniformRefine(); + + //=============================================// + // Setting loads and boundary conditions // + //=============================================// + + // source function, rhs + gsConstantFunction<> f(0.,0.,0.,3); + // surface load, neumann BC + gsConstantFunction<> g(0, 0, -1e8,3); + + // boundary conditions + gsBoundaryConditions<> bcInfo; + // Dirichlet BC are imposed separately for every component (coordinate) + for (index_t d = 0; d < 3; d++) + { + bcInfo.addCondition(0,boundary::west,condition_type::dirichlet,0,d); + } + // Neumann BC are imposed as one function + bcInfo.addCondition(1,boundary::east,condition_type::neumann,&g); + + //=============================================// + // Solving // + //=============================================// + + gsLinearMaterial materialMat(youngsModulus,poissonsRatio,3); + + // creating assembler + gsElasticityAssembler assembler(geometry,basis,bcInfo,f,&materialMat); + // gsElasticityAssembler assembler(geometry,basis,bcInfo,f); + assembler.options().setReal("YoungsModulus",youngsModulus); + assembler.options().setReal("PoissonsRatio",poissonsRatio); + assembler.options().setInt("MaterialLaw",materialLaw); + gsInfo << "Initialized system with " << assembler.numDofs() << " dofs.\n"; + + // setting Newton's method + gsIterative newton(assembler); + newton.options().setInt("MaxIters",50); + newton.options().setReal("AbsTol",1e-12); + newton.options().setInt("Verbosity",solver_verbosity::all); + newton.options().setInt("Solver",linear_solver::LDLT); + + gsInfo << "Solving...\n"; + gsStopwatch clock; + clock.restart(); + newton.solve(); + gsInfo << "Solved the system in " << clock.stop() <<"s.\n"; + + //=============================================// + // Output // + //=============================================// + + // solution to the nonlinear problem as an isogeometric displacement field + gsMultiPatch<> solutionNonlinear; + assembler.constructSolution(newton.solution(),newton.allFixedDofs(),solutionNonlinear); + // constructing stresses + gsPiecewiseFunction<> stresses; + assembler.constructCauchyStresses(solutionNonlinear,stresses,stress_components::von_mises); + + if (numPlotPoints > 0) + { + gsField<> nonlinearSolutionField(geometry,solutionNonlinear); + gsField<> stressField(assembler.patches(),stresses,true); + // creating a container to plot all fields to one Paraview file + std::map *> fields; + fields["Deformation"] = &nonlinearSolutionField; + fields["von Mises"] = &stressField; + gsWriteParaviewMultiPhysics(fields,"terrific_nle",numPlotPoints); + gsInfo << "Open \"terrific_nle.pvd\" in Paraview for visualization.\n"; + } + + return 0; +} diff --git a/examples/cooks_nonLinElast_2D.cpp b/examples/cooks_nonLinElast_2D.cpp index 6b0f884..277853b 100644 --- a/examples/cooks_nonLinElast_2D.cpp +++ b/examples/cooks_nonLinElast_2D.cpp @@ -7,6 +7,9 @@ #include #include #include +#include +#include +#include using namespace gismo; @@ -67,11 +70,16 @@ int main(int argc, char* argv[]){ // Solving // //=============================================// + gsNeoHookeLogMaterial materialMat(youngsModulus,poissonsRatio,2); + // gsLinearMaterial materialMat(youngsModulus,poissonsRatio); + // creating assembler - gsElasticityAssembler assembler(geometry,basisDisplacement,bcInfo,g); + gsElasticityAssembler assembler(geometry,basisDisplacement,bcInfo,g,&materialMat); + // gsElasticityAssembler assembler(geometry,basisDisplacement,bcInfo,g); assembler.options().setReal("YoungsModulus",youngsModulus); assembler.options().setReal("PoissonsRatio",poissonsRatio); assembler.options().setInt("MaterialLaw",material_law::neo_hooke_ln); + // assembler.options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff); gsInfo << "Initialized system with " << assembler.numDofs() << " dofs.\n"; // setting Newton's method diff --git a/examples/gsMaterial_test.cpp b/examples/gsMaterial_test.cpp new file mode 100644 index 0000000..68b6602 --- /dev/null +++ b/examples/gsMaterial_test.cpp @@ -0,0 +1,102 @@ +/** @file gsMaterialMatrix_test.cpp + + @brief Simple example for material matrix evaluations + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): H.M. Verhelst +*/ + +#include +#include +#include +#include + +using namespace gismo; + +int main (int argc, char** argv) +{ + gsCmdLine cmd("."); + + try { cmd.getValues(argc,argv); } catch (int rv) { return rv; } + + // Make geometry + gsMultiPatch<> mp, mp_def; + mp.addPatch( gsNurbsCreator<>::BSplineCube(1) ); // degree + mp.addAutoBoundaries(); + mp.uniformRefine(); // refines both + mp_def = mp; + mp_def.patch(0).coefs().col(0) *= 2; // deformation *2 in the x direction (col(0)) + + real_t E_modulus = 210e9; + real_t PoissonRatio = 0.3; + gsFunctionExpr<> E(std::to_string(E_modulus),3); + gsFunctionExpr<> nu(std::to_string(PoissonRatio),3); + + gsMatrix<> pts(3,6); + pts.col(0)<<0.125,0.375,0.5; + pts.col(1)<<0.375,0.125,0.5; + pts.col(2)<<0.125,0.25,0.5; + pts.col(3)<<0.25,0.125,0.5; + pts.col(4)<<0.25,0.25,0.5; + pts.col(5)<<0.5,0.5,0.5; + + gsMatrix<> Fres, Eres, Sres, Cres; + + gsLinearMaterial SvK(E_modulus, PoissonRatio, 3); + gsStopwatch time; + // ============================================================ + time.restart(); + // Using precomputation of material matrices + // Precompute the material matrix + gsMaterialData data; + SvK.precompute(&mp,&mp_def,0,pts,data); + SvK.eval_deformation_gradient_into(data,Fres); + SvK.eval_strain_into(data,Eres); + SvK.eval_stress_into(data,Sres); + SvK.eval_matrix_into(data,Cres); + gsInfo<<"Computation took "< SvK_F(&SvK, mp, mp_def); + gsMaterialEval SvK_E(&SvK, mp, mp_def); + gsMaterialEval SvK_S(&SvK, mp, mp_def); + gsMaterialEval SvK_C(&SvK, mp, mp_def); + SvK_F.piece(0).eval_into(pts,Fres); + SvK_E.piece(0).eval_into(pts,Eres); + SvK_S.piece(0).eval_into(pts,Sres); + SvK_C.piece(0).eval_into(pts,Cres); + gsInfo<<"Computation took "< #include #include +#include +#include using namespace gismo; @@ -75,8 +77,11 @@ int main(int argc, char* argv[]){ // Assembling & solving // //=============================================// + gsLinearMaterial materialMat(youngsModulus,poissonsRatio,2); + // creating assembler - gsElasticityAssembler assembler(geometry,basis,bcInfo,g); + // gsElasticityAssembler assembler(geometry,basis,bcInfo,g);//,materialMat); + gsElasticityAssembler assembler(geometry,basis,bcInfo,g,&materialMat); assembler.options().setReal("YoungsModulus",youngsModulus); assembler.options().setReal("PoissonsRatio",poissonsRatio); gsInfo<<"Assembling...\n"; diff --git a/examples/terrific_linElast_3D.cpp b/examples/terrific_linElast_3D.cpp index de6d769..f53d1e1 100644 --- a/examples/terrific_linElast_3D.cpp +++ b/examples/terrific_linElast_3D.cpp @@ -6,7 +6,8 @@ #include #include #include -#include +#include +#include using namespace gismo; @@ -72,8 +73,12 @@ int main(int argc, char* argv[]) // Assembling & solving // //=============================================// + gsLinearMaterial materialMat(youngsModulus,poissonsRatio,3); + // creating assembler - gsElasticityAssembler assembler(geometry,basis,bcInfo,f); + // gsElasticityAssembler assembler(geometry,basis,bcInfo,g);//,materialMat); + gsElasticityAssembler assembler(geometry,basis,bcInfo,f,&materialMat); + assembler.options().setReal("YoungsModulus",youngsModulus); assembler.options().setReal("PoissonsRatio",poissonsRatio); assembler.options().setInt("DirichletValues",dirichlet::l2Projection); diff --git a/examples/terrific_nonLinElast_3D.cpp b/examples/terrific_nonLinElast_3D.cpp index feb6cdf..fb92d06 100644 --- a/examples/terrific_nonLinElast_3D.cpp +++ b/examples/terrific_nonLinElast_3D.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include using namespace gismo; @@ -74,6 +75,8 @@ int main(int argc, char* argv[]){ // Solving // //=============================================// + gsLinearMaterial materialMat(youngsModulus,poissonsRatio,3); + // creating assembler gsElasticityAssembler assembler(geometry,basis,bcInfo,f); assembler.options().setReal("YoungsModulus",youngsModulus); diff --git a/gsALE.h b/src/gsALE.h similarity index 100% rename from gsALE.h rename to src/gsALE.h diff --git a/gsALE.hpp b/src/gsALE.hpp similarity index 99% rename from gsALE.hpp rename to src/gsALE.hpp index b9a65f4..3c578a5 100644 --- a/gsALE.hpp +++ b/src/gsALE.hpp @@ -42,7 +42,7 @@ gsALE::gsALE(gsMultiPatch & geometry, const gsMultiPatch & displacement for (gsMultiPatch<>::const_biterator it = geometry.bBegin(); it != geometry.bEnd(); ++it) for (index_t d = 0; d < geometry.parDim(); ++d) bcInfo.addCondition(it->patch,it->side(),condition_type::dirichlet,0,d); - gsConstantFunction rhs(gsVector<>::Zero(geometry.parDim()),geometry.parDim()); + gsConstantFunction rhs(gsVector::Zero(geometry.parDim()),geometry.parDim()); // define assembler according to the method if (methodALE == ale_method::TINE || methodALE == ale_method::TINE_StVK || methodALE == ale_method::ILE || methodALE == ale_method::LE) diff --git a/gsALE_.cpp b/src/gsALE_.cpp similarity index 100% rename from gsALE_.cpp rename to src/gsALE_.cpp diff --git a/gsBaseAssembler.h b/src/gsBaseAssembler.h similarity index 98% rename from gsBaseAssembler.h rename to src/gsBaseAssembler.h index 1d3d895..a2230c4 100644 --- a/gsBaseAssembler.h +++ b/src/gsBaseAssembler.h @@ -92,7 +92,7 @@ class gsBaseAssembler : public gsAssembler virtual void setRHS(const gsMatrix & rhs) {m_system.rhs() = rhs;} - virtual void setMatrix(const gsSparseMatrix & matrix) {m_system.matrix() = matrix;} + virtual void setMatrix(const gsSparseMatrix & mat) {m_system.matrix() = mat;} protected: using gsAssembler::m_pde_ptr; diff --git a/gsBaseAssembler.hpp b/src/gsBaseAssembler.hpp similarity index 100% rename from gsBaseAssembler.hpp rename to src/gsBaseAssembler.hpp diff --git a/gsBaseAssembler_.cpp b/src/gsBaseAssembler_.cpp similarity index 100% rename from gsBaseAssembler_.cpp rename to src/gsBaseAssembler_.cpp diff --git a/gsBasePde.h b/src/gsBasePde.h similarity index 100% rename from gsBasePde.h rename to src/gsBasePde.h diff --git a/gsBaseUtils.h b/src/gsBaseUtils.h similarity index 100% rename from gsBaseUtils.h rename to src/gsBaseUtils.h diff --git a/gsBiharmonicAssembler.h b/src/gsBiharmonicAssembler.h similarity index 100% rename from gsBiharmonicAssembler.h rename to src/gsBiharmonicAssembler.h diff --git a/gsBiharmonicAssembler.hpp b/src/gsBiharmonicAssembler.hpp similarity index 100% rename from gsBiharmonicAssembler.hpp rename to src/gsBiharmonicAssembler.hpp diff --git a/gsBiharmonicAssembler_.cpp b/src/gsBiharmonicAssembler_.cpp similarity index 100% rename from gsBiharmonicAssembler_.cpp rename to src/gsBiharmonicAssembler_.cpp diff --git a/gsElPoissonAssembler.h b/src/gsElPoissonAssembler.h similarity index 100% rename from gsElPoissonAssembler.h rename to src/gsElPoissonAssembler.h diff --git a/gsElPoissonAssembler.hpp b/src/gsElPoissonAssembler.hpp similarity index 100% rename from gsElPoissonAssembler.hpp rename to src/gsElPoissonAssembler.hpp diff --git a/gsElPoissonAssembler_.cpp b/src/gsElPoissonAssembler_.cpp similarity index 100% rename from gsElPoissonAssembler_.cpp rename to src/gsElPoissonAssembler_.cpp diff --git a/gsElTimeIntegrator.h b/src/gsElTimeIntegrator.h similarity index 100% rename from gsElTimeIntegrator.h rename to src/gsElTimeIntegrator.h diff --git a/gsElTimeIntegrator.hpp b/src/gsElTimeIntegrator.hpp similarity index 100% rename from gsElTimeIntegrator.hpp rename to src/gsElTimeIntegrator.hpp diff --git a/gsElTimeIntegrator_.cpp b/src/gsElTimeIntegrator_.cpp similarity index 100% rename from gsElTimeIntegrator_.cpp rename to src/gsElTimeIntegrator_.cpp diff --git a/gsElasticityAssembler.h b/src/gsElasticityAssembler.h similarity index 85% rename from gsElasticityAssembler.h rename to src/gsElasticityAssembler.h index d33356a..8bfd301 100644 --- a/gsElasticityAssembler.h +++ b/src/gsElasticityAssembler.h @@ -18,6 +18,8 @@ #include #include #include +#include +#include namespace gismo { @@ -37,7 +39,7 @@ class gsElasticityAssembler : public gsBaseAssembler gsElasticityAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, const gsBoundaryConditions & bconditions, - const gsFunction & body_force); + const gsFunction & body_force); /// @brief Constructor of mixed formulation (displacement + pressure) gsElasticityAssembler(const gsMultiPatch & patches, @@ -46,6 +48,20 @@ class gsElasticityAssembler : public gsBaseAssembler const gsBoundaryConditions & bconditions, const gsFunction & body_force); + /// @brief Constructor for displacement formulation + gsElasticityAssembler(const gsMultiPatch & patches, + const gsMultiBasis & basis, + const gsBoundaryConditions & bconditions, + const gsFunction & body_force, + const gsMaterialBase * material); + + /// @brief Constructor for displacement formulation + gsElasticityAssembler(const gsMultiPatch & patches, + const gsMultiBasis & basis, + const gsBoundaryConditions & bconditions, + const gsFunction & body_force, + gsMaterialContainer materials); + /// @brief Returns the list of default options for assembly static gsOptionList defaultOptions(); @@ -124,12 +140,17 @@ class gsElasticityAssembler : public gsBaseAssembler /// parametric dim = physical dim = deformation dim short_t m_dim; + // const std::vector::uPtr > m_pars; + using Base::m_pde_ptr; using Base::m_bases; using Base::m_ddof; using Base::m_options; using Base::m_system; using Base::eliminationMatrix; + + gsMaterialContainer m_materials; + }; diff --git a/gsElasticityAssembler.hpp b/src/gsElasticityAssembler.hpp similarity index 84% rename from gsElasticityAssembler.hpp rename to src/gsElasticityAssembler.hpp index 6e24899..1394c40 100644 --- a/gsElasticityAssembler.hpp +++ b/src/gsElasticityAssembler.hpp @@ -24,9 +24,11 @@ // Element visitors #include +#include #include #include #include +#include #include namespace gismo @@ -37,43 +39,67 @@ gsElasticityAssembler::gsElasticityAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, const gsBoundaryConditions & bconditions, const gsFunction & body_force) +: +gsElasticityAssembler(patches,basis,bconditions,body_force,gsMaterialContainer()) +{} + +template +gsElasticityAssembler::gsElasticityAssembler(const gsMultiPatch & patches, + gsMultiBasis const & basisDisp, + gsMultiBasis const & basisPres, + const gsBoundaryConditions & bconditions, + const gsFunction & body_force) +: +m_materials((index_t)0) { - // Originally concieved as a meaningful class, now gsPde is just a container for - // the domain, boundary conditions and the right-hand side; - // TUDO: change/remove gsPde from gsAssembler logic + // same as above gsPiecewiseFunction rightHandSides; rightHandSides.addPiece(body_force); typename gsPde::Ptr pde( new gsBasePde(patches,bconditions,rightHandSides) ); - // gsAssembler<>::initialize requires a vector of bases, one for each unknown; - // different bases are used to compute Dirichlet DoFs; - // but always the first basis is used for the assembly; - // TODO: change gsAssembler logic + // same as above m_dim = body_force.targetDim(); for (short_t d = 0; d < m_dim; ++d) - m_bases.push_back(basis); + m_bases.push_back(basisDisp); + m_bases.push_back(basisPres); Base::initialize(pde, m_bases, defaultOptions()); + m_options.setInt("MaterialLaw",material_law::mixed_hooke); } template -gsElasticityAssembler::gsElasticityAssembler(gsMultiPatch const & patches, - gsMultiBasis const & basisDisp, - gsMultiBasis const & basisPres, - gsBoundaryConditions const & bconditions, - const gsFunction & body_force) +gsElasticityAssembler::gsElasticityAssembler(const gsMultiPatch & patches, + const gsMultiBasis & basis, + const gsBoundaryConditions & bconditions, + const gsFunction & body_force, + const gsMaterialBase * material) +: +gsElasticityAssembler(patches,basis,bconditions,body_force,gsMaterialContainer(material,patches.nPatches())) +{} + +template +gsElasticityAssembler::gsElasticityAssembler(const gsMultiPatch & patches, + const gsMultiBasis & basis, + const gsBoundaryConditions & bconditions, + const gsFunction & body_force, + gsMaterialContainer materials) +: +m_materials(materials) { - // same as above + // Originally concieved as a meaningful class, now gsPde is just a container for + // the domain, boundary conditions and the right-hand side; + // TUDO: change/remove gsPde from gsAssembler logic gsPiecewiseFunction rightHandSides; rightHandSides.addPiece(body_force); typename gsPde::Ptr pde( new gsBasePde(patches,bconditions,rightHandSides) ); - // same as above + // gsAssembler<>::initialize requires a vector of bases, one for each unknown; + // different bases are used to compute Dirichlet DoFs; + // but always the first basis is used for the assembly; + // TODO: change gsAssembler logic m_dim = body_force.targetDim(); for (short_t d = 0; d < m_dim; ++d) - m_bases.push_back(basisDisp); - m_bases.push_back(basisPres); + m_bases.push_back(basis); Base::initialize(pde, m_bases, defaultOptions()); - m_options.setInt("MaterialLaw",material_law::mixed_hooke); } template @@ -147,7 +173,7 @@ void gsElasticityAssembler::assemble(bool saveEliminationMatrix) if (m_bases.size() == unsigned(m_dim)) // displacement formulation { GISMO_ENSURE(m_options.getInt("MaterialLaw") == material_law::hooke || m_options.getInt("MaterialLaw") == material_law::saint_venant_kirchhoff, - "Material law not specified OR not supported!"); + "Material law not specified OR not supported! Material law = "<::assemble(bool saveEliminationMatrix) eliminationMatrix.reservePerColumn(m_system.numColNz(m_bases[0],m_options)); } - gsVisitorLinearElasticity visitor(*m_pde_ptr, saveEliminationMatrix ? &eliminationMatrix : nullptr); - Base::template push >(visitor); + // if !composite + + // else + + + if (m_materials.size()!=0) + { + gsVisitorLinearElasticityMM visitor(*m_pde_ptr,m_materials, saveEliminationMatrix ? &eliminationMatrix : nullptr); + Base::template push >(visitor); + } + else + { + gsVisitorLinearElasticity visitor(*m_pde_ptr, saveEliminationMatrix ? &eliminationMatrix : nullptr); + Base::template push >(visitor); + } if (saveEliminationMatrix) { @@ -205,15 +244,27 @@ void gsElasticityAssembler::assemble(const gsMultiPatch & displacement) { GISMO_ENSURE(m_options.getInt("MaterialLaw") == material_law::saint_venant_kirchhoff || m_options.getInt("MaterialLaw") == material_law::neo_hooke_ln || - m_options.getInt("MaterialLaw") == material_law::neo_hooke_quad, + m_options.getInt("MaterialLaw") == material_law::neo_hooke_quad || + m_materials.size()!=0, "Material law not specified OR not supported!"); m_system.matrix().setZero(); reserve(); m_system.rhs().setZero(); // Compute volumetric integrals and write to the global linear system - gsVisitorNonLinearElasticity visitor(*m_pde_ptr,displacement); - Base::template push >(visitor); + if (m_materials.size()!=0) + { + // gsVisitorNonLinearElasticity visitor(*m_pde_ptr,displacement); + // Base::template push >(visitor); + gsVisitorNonLinearElasticityMM visitor(*m_pde_ptr,m_materials,displacement); + Base::template push >(visitor); + } + else + { + gsVisitorNonLinearElasticity visitor(*m_pde_ptr,displacement); + Base::template push >(visitor); + } + // Compute surface integrals and write to the global rhs vector // change to reuse rhs from linear system Base::template push >(m_pde_ptr->bc().neumannSides()); diff --git a/gsElasticityAssembler_.cpp b/src/gsElasticityAssembler_.cpp similarity index 100% rename from gsElasticityAssembler_.cpp rename to src/gsElasticityAssembler_.cpp diff --git a/gsElasticityFunctions.h b/src/gsElasticityFunctions.h similarity index 100% rename from gsElasticityFunctions.h rename to src/gsElasticityFunctions.h diff --git a/gsElasticityFunctions.hpp b/src/gsElasticityFunctions.hpp similarity index 100% rename from gsElasticityFunctions.hpp rename to src/gsElasticityFunctions.hpp diff --git a/gsElasticityFunctions_.cpp b/src/gsElasticityFunctions_.cpp similarity index 100% rename from gsElasticityFunctions_.cpp rename to src/gsElasticityFunctions_.cpp diff --git a/gsGeoUtils.h b/src/gsGeoUtils.h similarity index 100% rename from gsGeoUtils.h rename to src/gsGeoUtils.h diff --git a/gsGeoUtils.hpp b/src/gsGeoUtils.hpp similarity index 100% rename from gsGeoUtils.hpp rename to src/gsGeoUtils.hpp diff --git a/gsGeoUtils_.cpp b/src/gsGeoUtils_.cpp similarity index 100% rename from gsGeoUtils_.cpp rename to src/gsGeoUtils_.cpp diff --git a/gsIterative.h b/src/gsIterative.h similarity index 100% rename from gsIterative.h rename to src/gsIterative.h diff --git a/gsIterative.hpp b/src/gsIterative.hpp similarity index 100% rename from gsIterative.hpp rename to src/gsIterative.hpp diff --git a/gsIterative_.cpp b/src/gsIterative_.cpp similarity index 100% rename from gsIterative_.cpp rename to src/gsIterative_.cpp diff --git a/src/gsLinearMaterial.h b/src/gsLinearMaterial.h new file mode 100644 index 0000000..045ddef --- /dev/null +++ b/src/gsLinearMaterial.h @@ -0,0 +1,133 @@ +/** @file gsLinearMaterial.h + + @brief Provides a linear elastic material model + @todo check equation: + \f{align*}{ + \Psi(\mathbf{F}) &= \frac{\lambda}{2} \text{tr}^2(\mathbf{F}) + \mu \text{tr}(\mathbf{F}^T \mathbf{F})\\ + \mathbf{S} &= \lambda \text{tr}(\mathbf{E}) \mathbf{I} + 2\mu \mathbf{E}\\ + \mathbf{C} &= \lambda \mathbf{C}^{\text{vol}} + \mu \mathbf{C}^{\text{dev}} + \f} + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + H.M.Verhelst (2019 - ...., TU Delft) +*/ + +#pragma once + +#include +#include +#include + +namespace gismo +{ + +/** + * @brief The gsLinearMaterial class provides a linear elastic material model + * @ingroup Elasticity + * @tparam T + */ +template +class gsLinearMaterial : public gsMaterialBase +{ + +public: + using Base = gsMaterialBase; + + /** + * @brief Constructor with constant parameters + * @param E Young's modulus + * @param nu Poisson's ratio + * @param dim Dimension of the problem + */ + gsLinearMaterial( const T E, + const T nu, + const size_t dim) + : + gsLinearMaterial(gsConstantFunction(E,dim),gsConstantFunction(nu,dim)) + { + } + + /** + * @brief Constructor with function parameters + * @param E Young's modulus + * @param nu Poisson's ratio + */ + gsLinearMaterial(const gsFunctionSet & E, + const gsFunctionSet & nu) + : + Base() + { + this->setParameter(0,E); + this->setParameter(1,nu); + } + + /// See \ref gsMaterialBase.h for more details + void eval_stress_into(const gsMaterialData & data, gsMatrix & Sresult) const + { + const short_t dim = data.dim; + const index_t N = data.size; + + // Resize the result + Sresult.resize(dim*dim,N); + + // Lamé parameters + T E, nu; + T lambda, mu; + gsMatrix I = gsMatrix::Identity(dim,dim); + + for (index_t i=0; i!=N; i++) + { + E = data.m_parmat(0,i); + nu= data.m_parmat(1,i); + lambda = E * nu / ( ( 1. + nu ) * ( 1. - 2. * nu ) ); + mu = E / ( 2. * ( 1. + nu ) ); + + gsAsMatrix Emat = data.m_E.reshapeCol(i,dim,dim); + gsAsMatrix S = Sresult.reshapeCol(i,dim,dim); + S = lambda*Emat.trace()*I + 2*mu*Emat; + } + } + + /// See \ref gsMaterialBase.h for more details + void eval_matrix_into(const gsMaterialData & data, gsMatrix & Cresult) const + { + const short_t dim = data.dim; + const index_t N = data.size; + + // Voigt-size of the tensor + const index_t sz = (dim+1)*dim/2; + + // Resize the result + Cresult.resize(sz*sz,N); + + // Identity tensor + gsMatrix I = gsMatrix::Identity(dim,dim); + + // C tensors + gsMatrix Clambda, Cmu; + matrixTraceTensor(Clambda,I,I); + symmetricIdentityTensor(Cmu,I); + + // Lamé parameters + T E, nu; + T lambda, mu; + for (index_t i=0; i!=N; i++) + { + E = data.m_parmat(0,i); + nu= data.m_parmat(1,i); + lambda = E * nu / ( ( 1. + nu ) * ( 1. - 2. * nu ) ); + mu = E / ( 2. * ( 1. + nu ) ); + // Compute C + Cresult.reshapeCol(i,sz,sz) = lambda*Clambda + mu*Cmu; + } + } + +}; + +} diff --git a/gsMassAssembler.h b/src/gsMassAssembler.h similarity index 100% rename from gsMassAssembler.h rename to src/gsMassAssembler.h diff --git a/gsMassAssembler.hpp b/src/gsMassAssembler.hpp similarity index 100% rename from gsMassAssembler.hpp rename to src/gsMassAssembler.hpp diff --git a/gsMassAssembler_.cpp b/src/gsMassAssembler_.cpp similarity index 100% rename from gsMassAssembler_.cpp rename to src/gsMassAssembler_.cpp diff --git a/src/gsMaterialBase.h b/src/gsMaterialBase.h new file mode 100644 index 0000000..1fca860 --- /dev/null +++ b/src/gsMaterialBase.h @@ -0,0 +1,698 @@ +/** @file gsMaterialBase.h + + @brief Provides the base class for material models + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + O. Weeger (2012 - 2015, TU Kaiserslautern), + A.Shamanskiy (2016 - 2020, TU Kaiserslautern), + H.M.Verhelst (2019 - ...., TU Delft) +*/ + +#pragma once + +#include +#include +#include +#include + +namespace gismo +{ + +/* + NOTE ON THREADING: + The gsMaterialEvalSingle and gsMaterialEval classes contain a thread-safe data object m_data + which is used to store the material data for each thread. This is done to avoid + + */ + +template +class gsMaterialData; + +/** + * @brief Base class for material models + * @tparam T + * @ingroup Elasticity + */ +template +class gsMaterialBase +{ +public: + + typedef typename gsFunctionSet::Ptr function_ptr; + + /// Shared pointer for gsMaterialBase + typedef memory::shared_ptr< gsMaterialBase > Ptr; + + /// Unique pointer for gsMaterialBase + typedef memory::unique_ptr< gsMaterialBase > uPtr; + + /// Default constructor + gsMaterialBase() + : + gsMaterialBase(nullptr) + { + } + + /** + * @brief Constructor with a given density function + * @param Density Density function + */ + gsMaterialBase( const gsFunctionSet & Density) + : + gsMaterialBase(memory::make_shared(Density.clone().release())) + { + } + + /** + * @brief Constructor with a given density function + * @param Density Density function + */ + gsMaterialBase( const gsFunctionSet * Density) + : + gsMaterialBase(memory::make_shared_not_owned(Density)) + { + } + + /** + * @brief Constructor with a given density function + * @param Density Density function + */ + gsMaterialBase( const function_ptr & Density) + : + m_density(Density) + { + } + + GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION(gsMaterialBase, clone) + + /// Destructor + virtual ~gsMaterialBase() {}; + + /// Copy constructor (makes deep copy) + gsMaterialBase( const gsMaterialBase & other ) + { + operator=(other); + } + + /// Move constructor + gsMaterialBase( gsMaterialBase && other ) + { + operator=(give(other)); + } + + /// Copy assignment + gsMaterialBase & operator= ( const gsMaterialBase & other ) + { + if (this!=&other) + { + m_options = other.m_options; + + m_pars = other.m_pars; + m_density = other.m_density; + } + return *this; + } + + /// Move assignment + gsMaterialBase & operator= ( gsMaterialBase && other ) + { + m_options = give(other.m_options); + + + m_pars = give(other.m_pars); + m_density = give(other.m_density); + return *this; + } + + /** + * @brief Sets the default options + */ + virtual inline void defaultOptions() + { } + + /** + * @brief Returns the options + * + * @return \ref gsOptionList with the class options + */ + virtual inline gsOptionList & options() { return m_options; } + + /** + * @brief Sets the options + * + * @param[in] opt \ref gsOptionList + */ + virtual void setOptions(gsOptionList opt) {m_options.update(opt,gsOptionList::addIfUnknown); } + + /** + * @brief Precomputes the material data. + * + * @param[in] undeformed The undeformed configuration + * @param[in] patch The patch index + * @param[in] u The displacement field + * @param[out] data The material data + */ + virtual void precompute(const gsFunctionSet * undeformed, + const index_t patch, + const gsMatrix & u, + gsMaterialData & data) const + { + gsMapData map_ori, map_def; + this->_computeGeometricData(undeformed,nullptr,patch,u,map_ori,map_def); + this->_computeStrains(map_ori,map_def,data); + } + + /** + * @brief Precomputes the material data. + * + * @param[in] undeformed The undeformed configuration + * @param[in] deformed The deformed configuration + * @param[in] patch The patch index + * @param[in] u The displacement field + * @param[out] data The material data + */ + virtual void precompute(const gsFunctionSet * undeformed, + const gsFunctionSet * deformed, + const index_t patch, + const gsMatrix & u, + gsMaterialData & data) const + { + gsMapData map_ori, map_def; + this->_computeGeometricData(undeformed,deformed,patch,u,map_ori,map_def); + this->_computeStrains(map_ori,map_def,data); + } + + /** + * @brief Precomputes the material data from map data. + * + * @param[in] map_ori The undeformed map data + * @param[in] map_def The deformed map data + * @param[out] data The material data + */ + virtual void precompute(const gsMapData & map_ori, + const gsMapData & map_def, + gsMaterialData & data) const + { + this->_computeStrains(map_ori,map_def,data); + } + + /** + * @brief Evaluates the deformation gradient at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual void eval_deformation_gradient_into(const gsFunctionSet & undeformed, + const gsFunctionSet & deformed, + const index_t patch, + const gsMatrix& u, + gsMatrix & Fresult) const + { + // Compute map and parameters + gsMaterialData data; + gsMapData map_ori, map_def; + this->_computeGeometricData(&undeformed,&deformed,patch,u,map_ori,map_def); + this->_computeStrains(map_ori,map_def,data); + this->eval_deformation_gradient_into(data,Fresult); + } + + virtual void eval_deformation_gradient_into(const gsMaterialData & data, gsMatrix & Fresult) const + { + Fresult = data.m_F; + } + + + /** + * @brief Evaluates the strain tensor at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual void eval_strain_into( const gsFunctionSet & undeformed, + const gsFunctionSet & deformed, + const index_t patch, + const gsMatrix& u, + gsMatrix & Eresult) const + { + // Compute map and parameters + gsMaterialData data; + gsMapData map_ori, map_def; + this->_computeGeometricData(&undeformed,&deformed,patch,u,map_ori,map_def); + this->_computeStrains(map_ori,map_def,data); + this->eval_strain_into(data,Eresult); + } + + virtual void eval_strain_into(const gsMaterialData & data, gsMatrix & Eresult) const + { + Eresult = data.m_E; + } + + + /** + * @brief Evaluates the stress tensor at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual void eval_stress_into( const gsFunctionSet & undeformed, + const gsFunctionSet & deformed, + const index_t patch, + const gsMatrix& u, + gsMatrix & Sresult) const + { + // Compute map and parameters + gsMaterialData data; + gsMapData map_ori, map_def; + this->_computeGeometricData(&undeformed,&deformed,patch,u,map_ori,map_def); + this->_computeStrains(map_ori,map_def,data); + this->eval_stress_into(data,Sresult); + } + + virtual void eval_stress_into(const gsMaterialData & data, gsMatrix & Sresult) const + { GISMO_NO_IMPLEMENTATION; } + + /** + * @brief Evaluates the material tensor at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual void eval_matrix_into( const gsFunctionSet & undeformed, + const gsFunctionSet & deformed, + const index_t patch, + const gsMatrix& u, + gsMatrix & Cresult) const + { + // Compute map and parameters + gsMaterialData data; + gsMapData map_ori, map_def; + this->_computeGeometricData(&undeformed,&deformed,patch,u,map_ori,map_def); + this->_computeStrains(map_ori,map_def,data); + this->eval_matrix_into(data,Cresult); + } + + virtual void eval_matrix_into(const gsMaterialData & data, gsMatrix & Cresult) const + { GISMO_NO_IMPLEMENTATION; } + + /** + * @brief Evaluates the undamaged positive material tensor at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual void eval_matrix_pos_into( const gsFunctionSet & undeformed, + const gsFunctionSet & deformed, + const index_t patch, + const gsMatrix& u, + gsMatrix & Cresult) const + { + // Compute map and parameters + gsMaterialData data; + gsMapData map_ori, map_def; + this->_computeGeometricData(&undeformed,&deformed,patch,u,map_ori,map_def); + this->_computeStrains(map_ori,map_def,data); + this->eval_matrix_pos_into(data,Cresult); + } + + virtual void eval_matrix_pos_into(const gsMaterialData & data, gsMatrix & Cresult) const + { GISMO_NO_IMPLEMENTATION; } + + /** + * @brief Evaluates the energy at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual void eval_energy_into( const gsFunctionSet & undeformed, + const gsFunctionSet & deformed, + const index_t patch, + const gsMatrix& u, + gsMatrix & Presult) const + { + // Compute map and parameters + gsMaterialData data; + gsMapData map_ori, map_def; + this->_computeGeometricData(&undeformed,&deformed,patch,u,map_ori,map_def); + this->_computeStrains(map_ori,map_def,data); + this->eval_energy_into(data,Presult); + } + + virtual void eval_energy_into(const gsMaterialData & data, gsMatrix & Presult) const + { GISMO_NO_IMPLEMENTATION; } + + /// Sets the density + virtual inline void setDensity(function_ptr Density) { m_density = Density; } + /// Sets the density + virtual inline void setDensity(const gsFunctionSet & Density) + { + function_ptr fun = memory::make_shared(Density.clone().release()); + m_density = fun; + } + /// Returns true if a density is assigned + virtual inline bool hasDensity() const { return m_density!=nullptr; } + + /// Gets the Density + virtual const function_ptr getDensity() const {return m_density;} + + /** + * @brief Sets the material parameters. + * + * @param[in] pars Function pointers for the parameters in a container + */ + virtual inline void setParameters(const std::vector> &pars) + { + m_pars = pars; + } + + /** + * @brief Sets the material parameters. + * + * @param[in] pars Function pointers for the parameters in a container + * @param[in] parametric If true, the parameters are considered to be parametric + */ + virtual inline void setParameters(const std::vector &pars, bool parametric = true) + { + m_pars.resize(pars.size()); + for (size_t k = 0; k!=pars.size(); k++) + m_pars[k] = std::make_pair(pars[k],true); + } + + /** + * @brief Sets the material parameters. + * + * @param[in] pars Function pointers for the parameters in a container + */ + virtual inline void setParameter(const index_t i, const function_ptr &par, bool parametric = true) + { + if ((index_t)m_pars.size() < i+1) + m_pars.resize(i+1); + + m_pars[i] = std::make_pair(par,parametric); + } + + /** + * @brief Sets the material parameters. + * + * @param[in] pars Function pointers for the parameters in a container + */ + virtual inline void setParameters(const std::vector *,bool>> &pars) + { + m_pars.resize(pars.size()); + for (size_t k = 0; k!=pars.size(); k++) + m_pars[k] = std::make_pair(memory::make_shared_not_owned(pars[k].first),pars[k].second); + } + + /** + * @brief Sets the material parameters. + * + * @param[in] pars Function pointers for the parameters in a container + * @param[in] parametric If true, the parameters are considered to be parametric + */ + virtual inline void setParameters(const std::vector *> &pars, bool parametric = true) + { + m_pars.resize(pars.size()); + for (size_t k = 0; k!=pars.size(); k++) + m_pars[k] = std::make_pair(memory::make_shared_not_owned(pars[k]),parametric); + } + + /** + * @brief Sets the material parameters. + * + * @param[in] pars Function pointers for the parameters in a container + */ + virtual inline void setParameter(const index_t i, const gsFunctionSet &par, bool parametric = true) + { + if ((index_t)m_pars.size() < i+1) + m_pars.resize(i+1); + + m_pars[i] = std::make_pair(memory::make_shared(par.clone().release()),parametric); + } + + /** + * @brief Gets parameter i + * + * @param[in] i The parameter index + * + * @return The parameter. + */ + virtual inline const function_ptr getParameter(const index_t i) const + { + GISMO_ASSERT(i < (index_t)m_pars.size(),"Parameter "< & map_ori, + const gsMapData & map_def, + gsMaterialData & data) const + { + if (map_def.flags==0) + _computeStrains_impl(map_ori,data); + else + _computeStrains_impl(map_ori,map_def,data); + } + + /** + * @brief Computes the strains from undeformed map data only + * + * @param[in] map_ori The undeformed map data + * @param[out] data The material data + */ + void _computeStrains_impl(const gsMapData & map_ori, + gsMaterialData & data) const + { + GISMO_ASSERT(map_ori.flags & NEED_VALUE,"The undeformed map data should contain the point values"); + GISMO_ASSERT(map_ori.flags & NEED_JACOBIAN,"The undeformed map data should contain the Jacobian matrix"); + data.dim = map_ori.points.rows(); + data.size = map_ori.points.cols(); + data.patch = (map_ori.patchId==-1) ? 0 : map_ori.patchId; + data.m_F.resize(0,data.size); + data.m_E.resize(0,data.size); + data.m_parmat.resize(m_pars.size(),data.size); + + data.m_parmat.setZero(); + gsMatrix pars; + for (size_t v=0; v!=m_pars.size(); v++) + { + m_pars[v].first->piece(data.patch).eval_into((m_pars[v].second) ? + map_ori.points : + map_ori.values[0], + pars); + data.m_parmat.row(v) = pars; + } + } + + /** + * @brief Computes the strains from undeformed and deformed map data + * + * @param[in] map_ori The undeformed map data + * @param[in] map_def The deformed map data + * @param[out] data The material data + */ + void _computeStrains_impl(const gsMapData & map_ori, + const gsMapData & map_def, + gsMaterialData & data) const + { + GISMO_ASSERT(map_ori.flags & NEED_VALUE,"The undeformed map data should contain the point values"); + GISMO_ASSERT(map_ori.flags & NEED_JACOBIAN,"The undeformed map data should contain the Jacobian matrix"); + GISMO_ASSERT(map_def.flags & NEED_JACOBIAN,"The deformed map data should contain the Jacobian matrix"); + GISMO_ASSERT(map_ori.points.cols() == map_def.points.cols(),"The number of points in the undeformed and deformed maps should be the same"); + GISMO_ASSERT(map_ori.patchId == map_def.patchId,"The undeformed and deformed maps should be defined on the same patch"); + data.dim = map_ori.points.rows(); + data.size = map_ori.points.cols(); + data.patch = (map_ori.patchId==-1) ? 0 : map_ori.patchId; + data.m_F.resize(data.dim*data.dim,data.size); + data.m_E.resize(data.dim*data.dim,data.size); + data.m_parmat.resize(m_pars.size(),data.size); + + data.m_parmat.setZero(); + gsMatrix pars; + gsMatrix jac_ori, jac_def; + gsMatrix I = gsMatrix::Identity(data.dim,data.dim); + for (index_t k=0; k!= data.size; k++) + { + jac_ori = map_ori.jacobian(k); + jac_def = map_def.jacobian(k); + // deformation gradient F = I + du/dx + gsAsMatrix F = data.m_F.reshapeCol(k,data.dim,data.dim); + F = jac_def*(jac_ori.cramerInverse());; + // strain tensor E = 0.5*(F'*F-I), a.k.a. full geometric strain tensor + gsAsMatrix Emat = data.m_E.reshapeCol(k,data.dim,data.dim); + Emat = 0.5 * (F.transpose() * F - I); + } + + for (size_t v=0; v!=m_pars.size(); v++) + { + m_pars[v].first->piece(data.patch).eval_into(map_ori.values[0], pars); + data.m_parmat.row(v) = pars; + } + } + + /** + * @brief Computes the geometric data + * + * @param[in] undeformed The undeformed configuration + * @param[in] deformed The deformed configuration + * @param[in] patch The patch index + * @param[in] u The displacement field + * @param[out] map_ori The undeformed map data + * @param[out] map_def The deformed map data + */ + void _computeGeometricData(const gsFunctionSet * undeformed, + const gsFunctionSet * deformed, + const index_t patch, + const gsMatrix& u, + gsMapData & map_ori, + gsMapData & map_def) const + { + if (deformed==nullptr || deformed->nPieces()==0) + _computeGeometricData_impl(undeformed,patch,u,map_ori,map_def); + else + _computeGeometricData_impl(undeformed,deformed,patch,u,map_ori,map_def); + } + + /** + * @brief Computes the geometric data from undeformed configuration only + * + * @param[in] undeformed The undeformed configuration + * @param[in] patch The patch index + * @param[in] u The displacement field + * @param[out] map_ori The undeformed map data + * @param[out] map_def The deformed map data + */ + void _computeGeometricData_impl(const gsFunctionSet * undeformed, + const index_t patch, + const gsMatrix& u, + gsMapData & map_ori, + gsMapData & map_def) const + { + map_ori.clear(); + map_def.clear(); + GISMO_ASSERT(undeformed->domainDim()==u.rows(),"Geometric dimension and the point dimension are not the same!\n"<domainDim()<<"!="<function(patch).computeMap(map_ori); + } + + /** + * @brief Computes the geometric data from undeformed and deformed configurations + * + * @param[in] undeformed The undeformed configuration + * @param[in] deformed The deformed configuration + * @param[in] patch The patch index + * @param[in] u The displacement field + * @param[out] map_ori The undeformed map data + * @param[out] map_def The deformed map data + */ + void _computeGeometricData_impl(const gsFunctionSet * undeformed, + const gsFunctionSet * deformed, + const index_t patch, + const gsMatrix& u, + gsMapData & map_ori, + gsMapData & map_def) const + { + map_ori.clear(); + map_def.clear(); + GISMO_ASSERT(undeformed->domainDim()==deformed->domainDim(),"Geometric dimension and the point dimension are not the same!\n"<domainDim()<<"!="<domainDim()); + GISMO_ASSERT(undeformed->domainDim()==u.rows(),"Geometric dimension and the point dimension are not the same!\n"<domainDim()<<"!="<function(patch).computeMap(map_ori); + deformed->function(patch).computeMap(map_def); + } + +protected: + + std::vector< std::pair > m_pars; + function_ptr m_density; + + gsOptionList m_options; + + // Geometric data point + // mutable util::gsThreaded< gsMaterialData > m_data; + +public: + +# define Eigen gsEigen + EIGEN_MAKE_ALIGNED_OPERATOR_NEW //must be present whenever the class contains fixed size matrices +# undef Eigen + +}; + +/** + * @brief Material data container + * This class contains deformation gradients and strains + * @tparam T Real type + * @ingroup Elasticity + */ +template +class gsMaterialData +{ + +public: + + void membersSetZero() + { + m_F.setZero(); + m_E.setZero(); + m_rhoMat.setZero(); + } + + typename gsMaterialBase::function_ptr m_undeformed; + typename gsMaterialBase::function_ptr m_deformed; + + mutable gsMatrix m_parmat; + mutable gsMatrix m_rhoMat; + // mutable gsMatrix m_jac_ori, m_jac_def; + mutable gsMatrix m_F, m_E; + + mutable short_t dim; + mutable index_t size; + mutable index_t patch; +}; + +} diff --git a/src/gsMaterialBase_.cpp b/src/gsMaterialBase_.cpp new file mode 100644 index 0000000..a3bd431 --- /dev/null +++ b/src/gsMaterialBase_.cpp @@ -0,0 +1,11 @@ +// #include + +// #include +// #include + +// namespace gismo +// { +// CLASS_TEMPLATE_INST gsMaterialBase; + +// } + diff --git a/src/gsMaterialContainer.h b/src/gsMaterialContainer.h new file mode 100644 index 0000000..f02d665 --- /dev/null +++ b/src/gsMaterialContainer.h @@ -0,0 +1,196 @@ +/** @file gsMaterialContainer.h + + @brief Provides a container for material matrices for thin shells + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + H.M. Verhelst (2019-..., TU Delft) + A. Mantzaflaris (2019-..., Inria) +*/ + +#pragma once + +#include + +namespace gismo +{ + + +/** + * @brief This class serves as the evaluator of materials, based on \ref gsMaterialBase + * + * @tparam T Real type + * @tparam out Output type (see \ref gsMaterialOutput) + * + * @ingroup Elasticity + */ +template +class gsMaterialContainer // change name to PtrContainer +{ +public: + typedef typename std::vector::Ptr> Container; + + typedef typename Container::iterator iterator; + typedef typename Container::const_iterator const_iterator; + + /// Shared pointer for gsMaterialContainer + typedef memory::shared_ptr< gsMaterialContainer > Ptr; + + /// Unique pointer for gsMaterialContainer + typedef memory::unique_ptr< gsMaterialContainer > uPtr; +public: + + /// Constructor + gsMaterialContainer( index_t size = 0 ) + { + m_container.resize(size); + // To do: initialize with null pointers + } + + /// Constructor + gsMaterialContainer( const gsMaterialBase * mat, index_t size = 1 ) + { + m_container.reserve(size); + for (index_t i=0; i!=size; i++) + this->add(mat); + } + + gsMaterialContainer(const gsMaterialContainer & other) + { + // for (index_t k=0; k!=other.m_container.size(); k++) + // add(memory::make_unique(other.m_container.at(k))); + m_container = give(other.m_container); + } + + ~gsMaterialContainer() + { + // freeAll(m_container); + } + + /// Add a material matrix by copying argument + void add(const gsMaterialBase & mat) + { + m_container.push_back( memory::make_shared(mat.clone().release()) ); + } + + ///\brief Add a material matrix from a gsMaterialBase::uPtr + void add(const gsMaterialBase * mat) + { + m_container.push_back( memory::make_shared_not_owned(mat) ); + } + + /// Set a material matrix by copying argument + void set(const index_t i, const gsMaterialBase & mat) + { + m_container[i] = memory::make_shared(mat.clone().release()); + } + + ///\brief Set a material matrix from a gsMaterialBase::uPtr + void set(const index_t i, const gsMaterialBase * mat) + { + m_container[i] = memory::make_shared_not_owned(mat); + } + + ///\brief Set a material matrix from a gsMaterialBase::uPtr + void set(const index_t i, const typename gsMaterialBase::Ptr mat) + { + m_container[i] = mat; + } + + gsMaterialBase * piece(const index_t k) const + { + return m_container.at(k).get(); + } + + index_t size() const {return m_container.size();} + + std::ostream &print(std::ostream &os) const + { + os << "Material container with "< +// class gsXml< gsMaterialContainer > +// { +// private: +// gsXml() { } +// typedef gsMaterialContainer Object; + +// public: +// GSXML_COMMON_FUNCTIONS(gsMaterialContainer); +// static std::string tag () { return "MaterialContainer"; } +// static std::string type () { return ""; } + +// GSXML_GET_POINTER(Object); + +// static void get_into(gsXmlNode * node,Object & obj) +// { +// const int size = atoi(node->first_attribute("size")->value()); + +// // Read material inventory +// int count = countByTag("Material", node); +// std::vector::Ptr> mat(count); +// for (gsXmlNode * child = node->first_node("Material"); child; child = +// child->next_sibling("Material")) +// { +// const int i = atoi(child->first_attribute("index")->value()); +// mat[i] = memory::make_shared(gsXml>::get(child)); +// } + +// obj = gsMaterialContainer(size); +// for (gsXmlNode * child = node->first_node("group"); child; +// child = child->next_sibling("group")) +// { +// const int mIndex = atoi(child->first_attribute("material")->value()); +// std::istringstream group_str; +// group_str.str( child->value() ); + +// for(int patch; ( gsGetInt(group_str,patch)); ) +// obj.set(patch,mat[mIndex]); +// } + +// } + +// static gsXmlNode * put (const Object & obj, +// gsXmlTree & data) +// { +// GISMO_ERROR("Writing gsMaterialContainer to Xml is not implemented"); +// // gsWarn<<"Writing gsMaterialContainer to Xml is not implemented\n"; +// // gsXmlNode * result; +// // return result; +// // return putMaterial< Object >( obj,data ); +// } +// }; + +// } // namespace internal + +} // namespace gismo diff --git a/src/gsMaterialEval.h b/src/gsMaterialEval.h new file mode 100644 index 0000000..8ecb217 --- /dev/null +++ b/src/gsMaterialEval.h @@ -0,0 +1,387 @@ +/** @file gsMaterialEval.h + + @brief Evaluator for `gsMaterialBase` objects + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + H.M.Verhelst +*/ + +#pragma once + +#include +#include +#include +#include + +namespace gismo +{ + +template +class gsMaterialEvalSingle; + +/** + * @brief This class serves as the evaluator of materials, based on \ref gsMaterialBase + * @tparam T Real type + * @tparam out Output type (see \ref MaterialOutput) + * @tparam voigt Voigt notation (true: Voigt, false: tensor) + * @ingroup Elasticity + */ +template +class gsMaterialEval : public gsFunctionSet +{ + typedef typename gsMaterialBase::function_ptr function_ptr; + +public: + + /** + * @brief Constructor + * @param materialMatrices Material matrices + * @param undeformed Undeformed geometry + */ + gsMaterialEval( const gsMaterialContainer & materialMatrices, + const gsFunctionSet & undeformed) + : + m_materials(materialMatrices), + m_undeformed(undeformed), + m_deformed(nullptr) + { + this->_makePieces(m_undeformed); + } + + /** + * @brief Constructor + * @param materialMatrix Material matrix + * @param undeformed Undeformed geometry + */ + gsMaterialEval( gsMaterialBase * materialMatrix, + const gsFunctionSet & undeformed) + : + m_materials(undeformed.nPieces()), + m_undeformed(undeformed), + m_deformed(nullptr) + { + for (index_t p = 0; p!=undeformed.nPieces(); ++p) + m_materials.set(p,materialMatrix); + this->_makePieces(m_undeformed); + } + + /** + * @brief Constructor + * @param materialMatrices Material matrices + * @param undeformed Undeformed geometry + * @param deformed Deformed geometry + */ + gsMaterialEval( const gsMaterialContainer & materialMatrices, + const gsFunctionSet & undeformed, + const gsFunctionSet & deformed) + : + m_materials(materialMatrices), + m_undeformed(undeformed), + m_deformed(deformed) + { + this->_makePieces(m_undeformed,m_deformed); + } + + /** + * @brief Constructor + * @param materialMatrix Material matrix + * @param undeformed Undeformed geometry + * @param deformed Deformed geometry + */ + gsMaterialEval( gsMaterialBase * materialMatrix, + const gsFunctionSet & undeformed, + const gsFunctionSet & deformed) + : + m_materials(deformed.nPieces()), + m_undeformed(memory::make_shared_not_owned(undeformed.clone().release())), + m_deformed(memory::make_shared_not_owned(deformed.clone().release())) + { + for (index_t p = 0; p!=deformed.nPieces(); ++p) + m_materials.set(p,materialMatrix); + this->_makePieces(m_undeformed,m_deformed); + } + + /// Destructor + ~gsMaterialEval() + { + freeAll(m_pieces); + } + + /// Implementation of domainDimension, see \ref gsFunctionSet + short_t domainDim() const override {return this->piece(0).domainDim();} + + /** + * @brief Target dimension + * + * For a scalar (e.g. density) the target dimension is 1, for a vector (e.g. stress tensor in Voight notation) the target dimension is 3 and for a matrix (e.g. the material matrix) the target dimension is 9, which can be reshaped to a 3x3 matrix. + * + * @return Returns the target dimension depending on the specified type (scalar, vector, matrix etc.) + */ + short_t targetDim() const override { return this->piece(0).targetDim(); } + + /// Implementation of piece, see \ref gsFunction + const gsFunction & piece(const index_t p) const override + { + return *m_pieces[p]; + } + + /// Implementation of nPieces(), see \ref gsFunctionSet + index_t nPieces() const override {return m_pieces.size();} + + /// Implementation of eval_into, see \ref gsFunction + void eval_into(const gsMatrix& u, gsMatrix& result) const override + { GISMO_NO_IMPLEMENTATION; } + +protected: + /// Makes function pieces + void _makePieces(function_ptr undeformed) + { + m_pieces.resize(undeformed->nPieces()); + for (size_t p = 0; p!=m_pieces.size(); ++p) + m_pieces.at(p) = new gsMaterialEvalSingle(p,m_materials.piece(p),undeformed,nullptr); + } + + /// Makes function pieces + void _makePieces(function_ptr undeformed, function_ptr deformed) + { + m_pieces.resize(deformed->nPieces()); + for (size_t p = 0; p!=m_pieces.size(); ++p) + m_pieces.at(p) = new gsMaterialEvalSingle(p,m_materials.piece(p),undeformed,deformed); + } + +protected: + gsMaterialContainer m_materials; + function_ptr m_undeformed; + function_ptr m_deformed; + mutable std::vector *> m_pieces; +}; //gsMaterialEval + +/** + * @brief This class serves as the evaluator of materials, based on \ref gsMaterialBase + * + * @tparam T Real type + * @tparam out Output type (see \ref MaterialOutput) + * @tparam voigt Voigt notation (true: Voigt, false: tensor) + * + * @ingroup Elasticity + */ +template +class gsMaterialEvalSingle : public gsFunction +{ + typedef typename gsMaterialBase::function_ptr function_ptr; + +public: + + /** + * @brief Constructor + * @param patch Patch index + * @param materialMatrix Material matrix + * @param undeformed Undeformed geometry + * @param deformed Deformed geometry + */ + gsMaterialEvalSingle( index_t patch, + gsMaterialBase * materialMatrix, + function_ptr undeformed, + function_ptr deformed) + : + m_pIndex(patch), + m_material(materialMatrix), + m_undeformed(undeformed), + m_deformed(deformed), + m_dim(m_undeformed->domainDim()) + { + } + + /// Implementation of domainDimension, see \ref gsFunction + short_t domainDim() const override {return m_dim;} + + /** + * @brief Target dimension + * + * @return Returns the target dimension depending on the specified type (scalar, vector, matrix etc.) + */ + short_t targetDim() const override {return targetDim_impl();} + + +private: + + /// Implementation of \ref targetDim for density, energy + /// @todo Implement density + template + typename std::enable_if<_out==gsMaterialOutput::Psi, short_t>::type targetDim_impl() const + { + return 1; + }; + + /// Implementation of \ref targetDim for strain, stress + template + typename std::enable_if<(_out==gsMaterialOutput::F || + _out==gsMaterialOutput::E || + _out==gsMaterialOutput::S) && !_voigt , short_t>::type targetDim_impl() const + { + return m_dim*m_dim; + }; + + /// Implementation of \ref targetDim for strain, stress + template + typename std::enable_if<(_out==gsMaterialOutput::F || + _out==gsMaterialOutput::E || + _out==gsMaterialOutput::S) && _voigt , short_t>::type targetDim_impl() const + { + return m_dim*(m_dim+1)/2; + }; + + /// Implementation of \ref targetDim for matrix C (size = ((d+1)*d/2) x ((d+1)*d/2)) -- in Voigt + template + typename std::enable_if<_out==gsMaterialOutput::C && _voigt, short_t>::type targetDim_impl() const + { + return math::pow((m_dim+1)*m_dim/2,2); + }; + + /// Implementation of \ref targetDim for matrix C (size = ((d+1)*d/2) x ((d+1)*d/2)) + template + typename std::enable_if<_out==gsMaterialOutput::C && !_voigt, short_t>::type targetDim_impl() const + { + GISMO_NO_IMPLEMENTATION; + }; + + /// Implementation of \ref targetDim for matrix C (size = ((d+1)*d/2) x ((d+1)*d/2)) -- in Voigt + template + typename std::enable_if<_out==gsMaterialOutput::C_pos && _voigt, short_t>::type targetDim_impl() const + { + return math::pow((m_dim+1)*m_dim/2,2); + }; + + /// Implementation of \ref targetDim for matrix C (size = ((d+1)*d/2) x ((d+1)*d/2)) + template + typename std::enable_if<_out==gsMaterialOutput::C_pos && !_voigt, short_t>::type targetDim_impl() const + { + GISMO_NO_IMPLEMENTATION; + }; + +protected: + /// Sets the patch index + void setPatch(index_t p) {m_pIndex = p; } + +public: + + /// Implementation of eval_into, see \ref gsFunction + void eval_into(const gsMatrix& u, gsMatrix& result) const override + { + this->eval_into_impl(u,result); + } + +private: + /// Specialisation of \ref eval_into for density (TODO), energy + template + typename std::enable_if<_out==gsMaterialOutput::Psi, void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + m_material->eval_energy_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + } + + /// Specialisation of \ref eval_into for the deformation gradient + template + typename std::enable_if<_out==gsMaterialOutput::F && !_voigt, void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + m_material->eval_deformation_gradient_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + } + + /// Specialisation of \ref eval_into for the deformation gradient + template + typename std::enable_if<_out==gsMaterialOutput::F && _voigt , void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + gsMatrix tmp; + m_material->eval_deformation_gradient_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + result.resize(m_dim*(m_dim+1)/2,u.cols()); + calculate_voigt_strain(tmp, m_dim, result); + } + + /// Specialisation of \ref eval_into for the strain tensor + template + typename std::enable_if<_out==gsMaterialOutput::E && !_voigt, void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + m_material->eval_strain_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + } + + /// Specialisation of \ref eval_into for the strain tensor + template + typename std::enable_if<_out==gsMaterialOutput::E && _voigt , void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + gsMatrix tmp; + m_material->eval_strain_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + result.resize(m_dim*(m_dim+1)/2,u.cols()); + calculate_voigt_strain(tmp, m_dim, result); + } + + /// Specialisation of \ref eval_into for the stress tensor + template + typename std::enable_if<_out==gsMaterialOutput::S && !_voigt, void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + m_material->eval_stress_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + } + + + /// Specialisation of \ref eval_into for the strain tensor + template + typename std::enable_if<_out==gsMaterialOutput::S && _voigt, void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + gsMatrix tmp; + m_material->eval_stress_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + result.resize(m_dim*(m_dim+1)/2,u.cols()); + calculate_voigt_stress(tmp, m_dim, result); + } + + /// Specialisation of \ref eval_into for the material tensor + template + typename std::enable_if<_out==gsMaterialOutput::C && _voigt , void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + m_material->eval_matrix_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + } + + template + typename std::enable_if<_out==gsMaterialOutput::C && !_voigt , void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + GISMO_NO_IMPLEMENTATION; + } + + /// Specialisation of \ref eval_into for the material tensor + template + typename std::enable_if<_out==gsMaterialOutput::C_pos && _voigt , void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + m_material->eval_matrix_pos_into(*m_undeformed,*m_deformed,m_pIndex,u,result); + } + + template + typename std::enable_if<_out==gsMaterialOutput::C_pos && !_voigt , void>::type eval_into_impl(const gsMatrix& u, gsMatrix& result) const + { + GISMO_NO_IMPLEMENTATION; + } + +public: + + +public: + + /// Precompute the geometric data + void precompute(const gsMatrix& u) + { + m_material->precomputeData(m_pIndex,u); + } + + /// Get the geometric data + +protected: + index_t m_pIndex; + gsMaterialBase * m_material; + function_ptr m_undeformed; + function_ptr m_deformed; + const short_t m_dim; +}; + +} // namespace diff --git a/src/gsMaterialEval.hpp b/src/gsMaterialEval.hpp new file mode 100644 index 0000000..ba69e5a --- /dev/null +++ b/src/gsMaterialEval.hpp @@ -0,0 +1,22 @@ +/** @file gsMaterialEval.hpp + + @brief Provides an evaluator for materials + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + H.M. Verhelst +*/ + +#pragma once + +#include + +namespace gismo +{ + +} // end namespace diff --git a/src/gsMaterialEval_.cpp b/src/gsMaterialEval_.cpp new file mode 100644 index 0000000..6b2f0fd --- /dev/null +++ b/src/gsMaterialEval_.cpp @@ -0,0 +1,17 @@ +// #include + +// #include +// // #include + +// namespace gismo +// { +// CLASS_TEMPLATE_INST gsMaterialEval; +// CLASS_TEMPLATE_INST gsMaterialEval; +// CLASS_TEMPLATE_INST gsMaterialEval; +// CLASS_TEMPLATE_INST gsMaterialEval; + +// CLASS_TEMPLATE_INST gsMaterialEvalSingle; +// CLASS_TEMPLATE_INST gsMaterialEvalSingle; +// CLASS_TEMPLATE_INST gsMaterialEvalSingle; +// CLASS_TEMPLATE_INST gsMaterialEvalSingle; +// } \ No newline at end of file diff --git a/src/gsMaterialHelper.h b/src/gsMaterialHelper.h new file mode 100644 index 0000000..fbb3b41 --- /dev/null +++ b/src/gsMaterialHelper.h @@ -0,0 +1,392 @@ +/** @file gsMaterialHelper.h + + @brief + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + O. Weeger (2012 - 2015, TU Kaiserslautern), + A.Shamanskiy (2016 - 2020, TU Kaiserslautern), + H.M.Verhelst (2019 - ...., TU Delft) +*/ + +#pragma once + +#include +#include +#include +#include + +namespace gismo +{ + +template +class gsMaterialData; + + +template +class gsMaterialHelper +{ +public: + + typedef typename gsFunctionSet::Ptr function_ptr; + + /// Shared pointer for gsMaterialHelper + typedef memory::shared_ptr< gsMaterialHelper > Ptr; + + /// Unique pointer for gsMaterialHelper + typedef memory::unique_ptr< gsMaterialHelper > uPtr; + + + gsMaterialHelper() + : + gsMaterialHelper(nullptr,nullptr) + { + } + + gsMaterialHelper( const gsFunctionSet * mp) + : + gsMaterialHelper(mp,nullptr) + { + } + + gsMaterialHelper( const gsFunctionSet & mp) + : + gsMaterialHelper(&mp,nullptr) + { + } + + gsMaterialHelper( const gsFunctionSet & mp, + const gsFunctionSet & mp_def) + : + gsMaterialHelper(memory::make_shared(mp.clone().release()), + memory::make_shared(mp_def.clone().release())) + { + } + + gsMaterialHelper( const gsFunctionSet * mp, + const gsFunctionSet * mp_def) + : + gsMaterialHelper(memory::make_shared_not_owned(mp), + memory::make_shared_not_owned(mp_def)) + { + } + + gsMaterialHelper( const function_ptr & mp, + const function_ptr & mp_def) + { + m_undeformed.mine() = mp; + m_deformed.mine() = mp_def; + } + + GISMO_UPTR_FUNCTION_NO_IMPLEMENTATION(gsMaterialHelper, clone) + + /// Destructor + virtual ~gsMaterialHelper() {}; + + /// Copy constructor (makes deep copy) + gsMaterialHelper( const gsMaterialHelper & other ) + { + operator=(other); + } + + /// Move constructor + gsMaterialHelper( gsMaterialHelper && other ) + { + operator=(give(other)); + } + + gsMaterialHelper & operator= ( const gsMaterialHelper & other ) + { + if (this!=&other) + { + m_undeformed.mine() = other.m_undeformed.mine(); + m_deformed.mine() = other.m_deformed.mine(); + m_options = other.m_options; + } + return *this; + } + + gsMaterialHelper & operator= ( gsMaterialHelper && other ) + { + m_undeformed.mine() = give(other.m_undeformed.mine()); + m_deformed.mine() = give(other.m_deformed.mine()); + m_options = give(other.m_options); + return *this; + } + + /** + * @brief Sets the default options + */ + virtual inline void defaultOptions() + { } + + /** + * @brief Returns the options + * + * @return \ref gsOptionList with the class options + */ + virtual inline gsOptionList & options() { return m_options; } + + /** + * @brief Returns the dimension + * + * @return The dimension + */ + virtual inline short_t dim() const + { + GISMO_ASSERT(m_undeformed.mine()!=nullptr,"Undeformed geometry is not set!"); + return m_undeformed.mine()->domainDim(); + } + + /** + * @brief Sets the options + * + * @param[in] opt \ref gsOptionList + */ + virtual inline void setOptions(gsOptionList opt) {m_options.update(opt,gsOptionList::addIfUnknown); } + + /** + * @brief Evaluates the deformation gradient at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual inline void eval_deformation_gradient_and_strain_into(const index_t patch, const gsMatrix& u, gsMatrix & Fresult, gsMatrix & Eresult) const + { + const short_t dim = u.rows(); + // Compute map and parameters + this->_computeGeometricData(patch,u); + + const index_t N = u.cols(); + gsMatrix I = gsMatrix::Identity(dim,dim); + gsMatrix physDefJac; + Fresult.resize(dim*dim,N); + Eresult.resize(dim*dim,N); + for (index_t i=0; i!=N; i++) + { + gsAsMatrix jacdef = m_data.mine().m_jac_def.reshapeCol(i,dim,dim); + gsAsMatrix jacori = m_data.mine().m_jac_ori.reshapeCol(i,dim,dim); + + physDefJac = jacdef*(jacori.cramerInverse()); + // deformation gradient F = I + du/dx + gsAsMatrix F = Fresult.reshapeCol(i,dim,dim); + F = physDefJac; + + // strain tensor E = 0.5*(F'*F-I), a.k.a. full geometric strain tensor + gsAsMatrix E = Eresult.reshapeCol(i,dim,dim); + E = 0.5 * (F.transpose() * F - I); + } + } + + /** + * @brief Evaluates the deformation gradient at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual inline void eval_deformation_gradient_into(const index_t patch, const gsMatrix& u, gsMatrix & result) const + { + const short_t dim = u.rows(); + // Compute map and parameters + this->_computeGeometricData(patch,u); + + const index_t N = u.cols(); + gsMatrix I = gsMatrix::Identity(dim,dim); + gsMatrix physDefJac; + result.resize(dim*dim,N); + for (index_t i=0; i!=N; i++) + { + gsAsMatrix jacdef = m_data.mine().m_jac_def.reshapeCol(i,dim,dim); + gsAsMatrix jacori = m_data.mine().m_jac_ori.reshapeCol(i,dim,dim); + + physDefJac = jacdef*(jacori.cramerInverse()); + // deformation gradient F = I + du/dx + gsAsMatrix F = result.reshapeCol(i,dim,dim); + + F = physDefJac; + + } + } + + + /** + * @brief Evaluates the strain tensor at a set of given points + * + * @param[in] patch The patch + * @param[in] u The in-plane shell coordinates to be eveluated on + * + * @return + */ + virtual inline void eval_strain_into(const index_t patch, const gsMatrix& u, gsMatrix & result) const + { + const short_t dim = u.rows(); + const index_t N = u.cols(); + gsMatrix Fres; + // NOTE: Could be faster to compute F here in the loop. + this->eval_deformation_gradient_into(patch,u,Fres); + gsMatrix I = gsMatrix::Identity(dim,dim); + result.resize(dim*dim,N); + + for (index_t i=0; i!=N; i++) + { + // deformation gradient F = I + du/dx + gsAsMatrix F = Fres.reshapeCol(i,dim,dim); + // Green-Lagrange strain, E = 0.5*(F'*F-I), a.k.a. full geometric strain tensor + gsAsMatrix E = result.reshapeCol(i,dim,dim); + + E = 0.5 * (F.transpose() * F - I); + // E = 0.5 * ((F-I).transpose() + (F-I)); // small strains + } + } + + /// Returns true if a density is assigned + virtual inline bool hasDensity() const { return m_density!=nullptr; } + + virtual void setUndeformed(const gsFunctionSet * undeformed) + { + function_ptr f_ptr = memory::make_shared_not_owned(undeformed); + m_undeformed.mine() = f_ptr; + } + virtual void setDeformed(const gsFunctionSet * deformed) + { + function_ptr f_ptr = memory::make_shared_not_owned(deformed); + m_deformed.mine() = f_ptr; + } + + virtual void setUndeformed(const function_ptr undeformed) + { + m_undeformed.mine() = undeformed; + } + virtual void setDeformed(const function_ptr deformed) + { + m_deformed.mine() = deformed; + } + + const function_ptr getUndeformed() const { return m_undeformed.mine(); } + const function_ptr getDeformed() const { return m_deformed.mine(); } + + virtual bool hasUndeformed() const + { return m_deformed.mine()!=nullptr; } + virtual bool hasDeformed() const { return m_deformed.mine()!=nullptr; } + + void precompute(index_t patch, const gsMatrix& u) + { + this->computeGeometricData(patch,u); + } + + const gsMaterialData & data() const{ return m_data.mine(); } + + + /** + * + */ + const gsMaterialData & getData() const{ return m_data.mine(); } + + void setData(const gsMaterialData & data) { m_data.mine() = data; } + +protected: + void _computeGeometricData(index_t patch, const gsMatrix& u) const + { + const short_t dim = u.rows(); + GISMO_ASSERT(m_undeformed.mine()!=nullptr,"Undeformed function is not set"); + GISMO_ASSERT(m_deformed.mine()!=nullptr,"Deformed function is not set"); + GISMO_ASSERT(m_undeformed.mine()->domainDim()==dim,"Geometric dimension and the point dimension are not the same!"); + GISMO_ASSERT(m_deformed.mine()->domainDim()==dim,"Geometric dimension and the point dimension are not the same!"); + m_data.mine().m_jac_ori.resize(dim*dim,u.cols()); + m_data.mine().m_jac_def.resize(dim*dim,u.cols()); + + gsMapData map_ori, map_def; + map_def.flags = map_ori.flags = NEED_JACOBIAN; + map_def.points = map_ori.points = u; + m_undeformed.mine()->function(patch).computeMap(map_ori); + m_deformed.mine()->function(patch).computeMap(map_def); + + for (index_t k=0; k!= u.cols(); k++) + { + m_data.mine().m_jac_ori.reshapeCol(k,dim,dim) = map_ori.jacobian(k); + m_data.mine().m_jac_def.reshapeCol(k,dim,dim) = map_def.jacobian(k); + } + } + + void _computeParameterData(index_t patch, const gsMatrix& u) const + { + const short_t dim = u.rows(); + GISMO_ASSERT(m_undeformed.mine()!=nullptr,"Undeformed function is not set"); + GISMO_ASSERT(m_undeformed.mine()->domainDim()==dim,"Geometric dimension and the point dimension are not the same!"); + m_data.mine().m_parmat.resize(m_pars.size(),u.cols()); + + gsMatrix tmp; + + gsMapData map; + map.flags = NEED_VALUE; + map.points = u; + static_cast&>(m_undeformed.mine()->piece(patch) ).computeMap(map); + + m_data.mine().m_parmat.resize(m_pars.size(),map.values[0].cols()); + m_data.mine().m_parmat.setZero(); + + for (size_t v=0; v!=m_pars.size(); v++) + { + m_pars[v]->piece(patch).eval_into(map.values[0], tmp); + m_data.mine().m_parmat.row(v) = tmp; + } + } + + +protected: + + void membersSetZero() + { + m_data.mine().membersSetZero(); + } + +protected: + + util::gsThreaded m_undeformed; + util::gsThreaded m_deformed; + + std::vector< function_ptr > m_pars; + function_ptr m_density; + + gsOptionList m_options; + + // Geometric data point + mutable util::gsThreaded< gsMaterialData > m_data; + +public: + +# define Eigen gsEigen + EIGEN_MAKE_ALIGNED_OPERATOR_NEW //must be present whenever the class contains fixed size matrices +# undef Eigen + +}; + +template +class gsMaterialData +{ + +public: + + void membersSetZero() + { + m_jac_ori.setZero(); + m_jac_def.setZero(); + m_rhoMat.setZero(); + } + + mutable gsMatrix m_parmat; + mutable gsMatrix m_rhoMat; + mutable gsMatrix m_jac_ori, m_jac_def; +}; + + +} diff --git a/src/gsMaterialUtils.h b/src/gsMaterialUtils.h new file mode 100644 index 0000000..0fe2ce3 --- /dev/null +++ b/src/gsMaterialUtils.h @@ -0,0 +1,72 @@ +/** @file gsMaterialEval.h + + @brief Utilities for gsMaterial + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + H.M.Verhelst +*/ + +#pragma once + +namespace gismo +{ + +enum class gsMaterialOutput : short_t +{ + /// @brief Deformation gradient + F = 00, + /// @brief Deformation gradient in Voigt notation + F_voigt = 01, + /// @brief Strain + E = 10, + /// @brief Strain in Voigt notation + E_voigt = 11, + /// @brief Stress + S = 20, + /// @brief Stress in Voigt notation + S_voigt = 21, + /// @brief Material matrix in Voigt notation + C = 30, + /// @brief Material matrix in Voigt notation undegraded + C_pos = 31, + /// @brief Energy + Psi = 40, + /// @brief Energy + Psi_pos = 41, +}; + +template +void calculate_voigt_stress(const gsMatrix& tmp, short_t d, gsMatrix& result) +{ + GISMO_ASSERT(tmp.cols() % d == 0, "Invalid size of stress matrix"); + result.resize(d*(d+1)/2, tmp.cols()/d); + for (index_t k = 0; k < tmp.cols()/d; k++) + { + gsMatrix S = tmp.reshapeCol(k, d, d); // in tensor notation + gsVector S_voigt; + voigtStress(S_voigt, S); // Convert to Voigt notation + result.col(k) = S_voigt; + } +} + +template +void calculate_voigt_strain(const gsMatrix& tmp, short_t d, gsMatrix& result) +{ + GISMO_ASSERT(tmp.cols() % d == 0, "Invalid size of strain matrix"); + result.resize(d*(d+1)/2, tmp.cols()/d); + for (index_t k = 0; k < tmp.cols()/d; k++) + { + gsMatrix E = tmp.reshapeCol(k, d, d); // in tensor notation + gsVector E_voigt; + voigtStrain(E_voigt, E); // Convert to Voigt notation + result.col(k) = E_voigt; + } +} + +} // namespace diff --git a/src/gsMaterial_.cpp b/src/gsMaterial_.cpp new file mode 100644 index 0000000..a4ed776 --- /dev/null +++ b/src/gsMaterial_.cpp @@ -0,0 +1,14 @@ +#include + +#include +#include +#include + +namespace gismo +{ + CLASS_TEMPLATE_INST gsLinearMaterial; + CLASS_TEMPLATE_INST gsNeoHookeLogMaterial; + CLASS_TEMPLATE_INST gsNeoHookeQuadMaterial; + +} + diff --git a/src/gsMooneyRivlinMaterial.h b/src/gsMooneyRivlinMaterial.h new file mode 100644 index 0000000..2bb4f3a --- /dev/null +++ b/src/gsMooneyRivlinMaterial.h @@ -0,0 +1,139 @@ +/** @file gsMooneyRivlinMaterial.h + + @brief + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + O. Weeger (2012 - 2015, TU Kaiserslautern), + A.Shamanskiy (2016 - 2020, TU Kaiserslautern), + H.M.Verhelst (2019 - ...., TU Delft) +*/ + +#pragma once + +#include +#include +#include + +namespace gismo +{ + +template +class gsMooneyRivlinMaterial : public gsMaterialBase +{ + +public: + using Base = gsMaterialBase; + + gsMooneyRivlinMaterial( const T c1, + const T c2, + short_t d) + : + gsMooneyRivlinMaterial( gsConstantFunction(c1,d), + gsConstantFunction(c2,d)) + { + } + + gsMooneyRivlinMaterial( const gsFunctionSet & c1, + const gsFunctionSet & c2) + : + Base() + { + this->setParameter(0,c1); + this->setParameter(1,c2); + } + + void eval_stress_into(const gsMaterialData & data, gsMatrix & Sresult) const + { + GISMO_NO_IMPLEMENTATION; + // const short_t dim = data.dim; + // const index_t N = data.size; + // // Compute the deformation gradient and strain + // gsMatrix Fres, Eres; + // Base::eval_deformation_gradient_and_strain_into(patch,u,Fres,Eres); + + // // Resize the result + // result.resize(dim*dim,u.cols()); + + // // Lamé parameters + // T E, nu; + // T lambda, mu; + // gsMatrix I = gsMatrix::Identity(dim,dim); + // gsMatrix RCG, RCGinv; + // T J; + // for (index_t i=0; i!=u.cols(); i++) + // { + // E = m_data.mine().m_parmat(0,i); + // nu= m_data.mine().m_parmat(1,i); + // lambda = E * nu / ( ( 1. + nu ) * ( 1. - 2. * nu ) ); + // mu = E / ( 2. * ( 1. + nu ) ); + + // gsAsMatrix F = Fres.reshapeCol(i,dim,dim); + // gsAsMatrix E = Eres.reshapeCol(i,dim,dim); + // gsAsMatrix S = result.reshapeCol(i,dim,dim); + // J = F.determinant(); + // RCG = F.transpose() * F; + // RCGinv = RCG.cramerInverse(); + // S = (lambda*log(J)-mu)*RCGinv + mu*I; + // } + } + + void eval_matrix_into(const gsMaterialData & data, gsMatrix & Cresult) const + { + GISMO_NO_IMPLEMENTATION; + // const short_t dim = data.dim; + // const index_t N = data.size; + // // Compute the deformation gradient and strain + // gsMatrix Fres, Eres; + // Base::eval_deformation_gradient_and_strain_into(patch,u,Fres,Eres); + + // // Voigt-size of the tensor + // const index_t sz = (dim+1)*dim/2; + + // // Resize the result + // result.resize(sz*sz,u.cols()); + + // // Identity tensor + // gsMatrix I = gsMatrix::Identity(dim,dim); + // gsMatrix C, Ctemp; + // gsMatrix RCG, RCGinv; + // T J; + + // // C tensors + // gsMatrix Clambda, Cmu; + // matrixTraceTensor(Clambda,I,I); + // symmetricIdentityTensor(Cmu,I); + + // // Lamé parameters + // T E, nu; + // T lambda, mu; + // for (index_t i=0; i!=u.cols(); i++) + // { + // E = m_data.mine().m_parmat(0,i); + // nu= m_data.mine().m_parmat(1,i); + // lambda = E * nu / ( ( 1. + nu ) * ( 1. - 2. * nu ) ); + // mu = E / ( 2. * ( 1. + nu ) ); + + // gsAsMatrix F = Fres.reshapeCol(i,dim,dim); + // gsAsMatrix E = Eres.reshapeCol(i,dim,dim); + + // J = F.determinant(); + // RCG = F.transpose() * F; + // RCGinv = RCG.cramerInverse(); + + // // Compute C + // matrixTraceTensor(C,RCGinv,RCGinv); + // C *= lambda; + // symmetricIdentityTensor(Ctemp,RCGinv); + // C += (mu-lambda*log(J))*Ctemp; + // result.reshapeCol(i,sz,sz) = C; + // } + } +}; + +} diff --git a/gsMuscleAssembler.h b/src/gsMuscleAssembler.h similarity index 100% rename from gsMuscleAssembler.h rename to src/gsMuscleAssembler.h diff --git a/gsMuscleAssembler.hpp b/src/gsMuscleAssembler.hpp similarity index 100% rename from gsMuscleAssembler.hpp rename to src/gsMuscleAssembler.hpp diff --git a/gsMuscleAssembler_.cpp b/src/gsMuscleAssembler_.cpp similarity index 100% rename from gsMuscleAssembler_.cpp rename to src/gsMuscleAssembler_.cpp diff --git a/src/gsNeoHookeLogMaterial.h b/src/gsNeoHookeLogMaterial.h new file mode 100644 index 0000000..5d60eeb --- /dev/null +++ b/src/gsNeoHookeLogMaterial.h @@ -0,0 +1,150 @@ +/** @file gsNeoHookeLogMaterial.h + + @brief Provides a neo-Hookean material model with logarithmic strain + @todo check equation: + \f{align*}{ + \Psi(\mathbf{F}) &= \frac{\lambda}{2} \log^2(J) - \mu \log(J) + \frac{\mu}{2} \text{tr}(\mathbf{C}^{\text{vol}})\\ + \mathbf{S} &= \lambda \log(J) \mathbf{C}^{-1} + \mu \mathbf{I}\\ + \mathbf{C} &= \lambda \mathbf{C}^{\text{vol}} + \mu \mathbf{C}^{\text{dev}} + \f} + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + H.M.Verhelst (2019 - ...., TU Delft) +*/ + +#pragma once + +#include +#include + +namespace gismo +{ + +/** + * @brief The gsNeoHookeLogMaterial class provides a neo-Hookean material model with logarithmic strain + * @ingroup Elasticity + * @tparam T + */ +template +class gsNeoHookeLogMaterial : public gsMaterialBase +{ + +public: + using Base = gsMaterialBase; + + /** + * @brief Constructor with constant parameters + * @param E Young's modulus + * @param nu Poisson's ratio + * @param dim Dimension of the problem + */ + gsNeoHookeLogMaterial( const T E, + const T nu, + short_t d) + : + gsNeoHookeLogMaterial(gsConstantFunction(E,d), + gsConstantFunction(nu,d)) + { + } + + /** + * @brief Constructor with function parameters + * @param E Young's modulus + * @param nu Poisson's ratio + */ + gsNeoHookeLogMaterial( const gsFunctionSet & E, + const gsFunctionSet & nu) + : + Base() + { + this->setParameter(0,E); + this->setParameter(1,nu); + } + + /// See \ref gsMaterialBase.h for more details + void eval_stress_into(const gsMaterialData & data, gsMatrix & Sresult) const override + { + const short_t dim = data.dim; + const index_t N = data.size; + + // Resize the result + Sresult.resize(dim*dim,N); + + // Lamé parameters + T E, nu; + T lambda, mu; + gsMatrix I = gsMatrix::Identity(dim,dim); + gsMatrix RCG, RCGinv; + T J; + for (index_t i=0; i!=N; i++) + { + E = data.m_parmat(0,i); + nu= data.m_parmat(1,i); + lambda = E * nu / ( ( 1. + nu ) * ( 1. - 2. * nu ) ); + mu = E / ( 2. * ( 1. + nu ) ); + + gsAsMatrix F = data.m_F.reshapeCol(i,dim,dim); + gsAsMatrix S = Sresult.reshapeCol(i,dim,dim); + J = F.determinant(); + RCG = F.transpose() * F; + RCGinv = RCG.cramerInverse(); + S = (lambda*log(J)-mu)*RCGinv + mu*I; + } + } + + /// See \ref gsMaterialBase.h for more details + void eval_matrix_into(const gsMaterialData & data, gsMatrix & Cresult) const override + { + const short_t dim = data.dim; + const index_t N = data.size; + + // Voigt-size of the tensor + const index_t sz = (dim+1)*dim/2; + + // Resize the result + Cresult.resize(sz*sz,N); + + // Identity tensor + gsMatrix I = gsMatrix::Identity(dim,dim); + gsMatrix C, Ctemp; + gsMatrix RCG, RCGinv; + T J; + + // C tensors + gsMatrix Clambda, Cmu; + matrixTraceTensor(Clambda,I,I); + symmetricIdentityTensor(Cmu,I); + + // Lamé parameters + T E, nu; + T lambda, mu; + for (index_t i=0; i!=N; i++) + { + E = data.m_parmat(0,i); + nu= data.m_parmat(1,i); + lambda = E * nu / ( ( 1. + nu ) * ( 1. - 2. * nu ) ); + mu = E / ( 2. * ( 1. + nu ) ); + + gsAsMatrix F = data.m_F.reshapeCol(i,dim,dim); + + J = F.determinant(); + RCG = F.transpose() * F; + RCGinv = RCG.cramerInverse(); + + // Compute C + matrixTraceTensor(C,RCGinv,RCGinv); + C *= lambda; + symmetricIdentityTensor(Ctemp,RCGinv); + C += (mu-lambda*log(J))*Ctemp; + Cresult.reshapeCol(i,sz,sz) = C; + } + } +}; + +} diff --git a/src/gsNeoHookeQuadMaterial.h b/src/gsNeoHookeQuadMaterial.h new file mode 100644 index 0000000..337a0a7 --- /dev/null +++ b/src/gsNeoHookeQuadMaterial.h @@ -0,0 +1,145 @@ +/** @file gsNeoHookeQuadMaterial.h + + @brief Provides a neo-Hookean material model + @todo equation + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + H.M.Verhelst (2019 - ...., TU Delft) +*/ + +#pragma once + +#include +#include + +namespace gismo +{ + +/** + * @brief The gsNeoHookeQuadMaterial class provides a neo-Hookean material model + * @ingroup Elasticity + * @tparam T + */ +template +class gsNeoHookeQuadMaterial : public gsMaterialBase +{ + +public: + using Base = gsMaterialBase; + + /** + * @brief Constructor with constant parameters + * @param E Young's modulus + * @param nu Poisson's ratio + * @param dim Dimension of the problem + */ + gsNeoHookeQuadMaterial( const T E, + const T nu, + short_t d) + : + gsNeoHookeQuadMaterial(gsConstantFunction(E,d), + gsConstantFunction(nu,d)) + { + } + + /** + * @brief Constructor with function parameters + * @param E Young's modulus + * @param nu Poisson's ratio + */ + gsNeoHookeQuadMaterial( const gsFunctionSet & E, + const gsFunctionSet & nu) + : + Base() + { + this->setParameter(0,E); + this->setParameter(1,nu); + } + + /// See \ref gsMaterialBase.h for more details + void eval_stress_into(const gsMaterialData & data, gsMatrix & Sresult) const + { + const short_t dim = data.dim; + const index_t N = data.size; + + // Resize the result + Sresult.resize(dim*dim,N); + + // Lamé parameters + T E, nu; + T lambda, mu; + gsMatrix I = gsMatrix::Identity(dim,dim); + gsMatrix RCG, RCGinv; + T J; + for (index_t i=0; i!=N; i++) + { + E = data.m_parmat(0,i); + nu= data.m_parmat(1,i); + lambda = E * nu / ( ( 1. + nu ) * ( 1. - 2. * nu ) ); + mu = E / ( 2. * ( 1. + nu ) ); + + gsAsMatrix F = data.m_F.reshapeCol(i,dim,dim); + gsAsMatrix S = Sresult.reshapeCol(i,dim,dim); + J = F.determinant(); + RCG = F.transpose() * F; + RCGinv = RCG.cramerInverse(); + S = (lambda*(J*J-1)/2-mu)*RCGinv + mu*I; + } + } + + /// See \ref gsMaterialBase.h for more details + void eval_matrix_into(const gsMaterialData & data, gsMatrix & Cresult) const + { + const short_t dim = data.dim; + const index_t N = data.size; + + // Voigt-size of the tensor + const index_t sz = (dim+1)*dim/2; + + // Resize the result + Cresult.resize(sz*sz,N); + + // Identity tensor + gsMatrix I = gsMatrix::Identity(dim,dim); + gsMatrix C, Ctemp; + gsMatrix RCG, RCGinv; + T J; + + // C tensors + gsMatrix Clambda, Cmu; + matrixTraceTensor(Clambda,I,I); + symmetricIdentityTensor(Cmu,I); + + // Lamé parameters + T E, nu; + T lambda, mu; + for (index_t i=0; i!=N; i++) + { + E = data.m_parmat(0,i); + nu= data.m_parmat(1,i); + lambda = E * nu / ( ( 1. + nu ) * ( 1. - 2. * nu ) ); + mu = E / ( 2. * ( 1. + nu ) ); + + gsAsMatrix F = data.m_F.reshapeCol(i,dim,dim); + + J = F.determinant(); + RCG = F.transpose() * F; + RCGinv = RCG.cramerInverse(); + + // Compute C + matrixTraceTensor(C,RCGinv,RCGinv); + C *= lambda*J*J; + symmetricIdentityTensor(Ctemp,RCGinv); + C += (mu-lambda*(J*J-1)/2)*Ctemp; + Cresult.reshapeCol(i,sz,sz) = C; + } + } +}; + +} diff --git a/gsNsAssembler.h b/src/gsNsAssembler.h similarity index 100% rename from gsNsAssembler.h rename to src/gsNsAssembler.h diff --git a/gsNsAssembler.hpp b/src/gsNsAssembler.hpp similarity index 100% rename from gsNsAssembler.hpp rename to src/gsNsAssembler.hpp diff --git a/gsNsAssembler_.cpp b/src/gsNsAssembler_.cpp similarity index 100% rename from gsNsAssembler_.cpp rename to src/gsNsAssembler_.cpp diff --git a/gsNsTimeIntegrator.h b/src/gsNsTimeIntegrator.h similarity index 100% rename from gsNsTimeIntegrator.h rename to src/gsNsTimeIntegrator.h diff --git a/gsNsTimeIntegrator.hpp b/src/gsNsTimeIntegrator.hpp similarity index 100% rename from gsNsTimeIntegrator.hpp rename to src/gsNsTimeIntegrator.hpp diff --git a/gsNsTimeIntegrator_.cpp b/src/gsNsTimeIntegrator_.cpp similarity index 100% rename from gsNsTimeIntegrator_.cpp rename to src/gsNsTimeIntegrator_.cpp diff --git a/gsPartitionedFSI.h b/src/gsPartitionedFSI.h similarity index 100% rename from gsPartitionedFSI.h rename to src/gsPartitionedFSI.h diff --git a/gsPartitionedFSI.hpp b/src/gsPartitionedFSI.hpp similarity index 100% rename from gsPartitionedFSI.hpp rename to src/gsPartitionedFSI.hpp diff --git a/gsPartitionedFSI_.cpp b/src/gsPartitionedFSI_.cpp similarity index 100% rename from gsPartitionedFSI_.cpp rename to src/gsPartitionedFSI_.cpp diff --git a/gsThermoAssembler.h b/src/gsThermoAssembler.h similarity index 100% rename from gsThermoAssembler.h rename to src/gsThermoAssembler.h diff --git a/gsThermoAssembler.hpp b/src/gsThermoAssembler.hpp similarity index 100% rename from gsThermoAssembler.hpp rename to src/gsThermoAssembler.hpp diff --git a/gsThermoAssembler_.cpp b/src/gsThermoAssembler_.cpp similarity index 100% rename from gsThermoAssembler_.cpp rename to src/gsThermoAssembler_.cpp diff --git a/gsVisitorBiharmonic.h b/src/gsVisitorBiharmonic.h similarity index 100% rename from gsVisitorBiharmonic.h rename to src/gsVisitorBiharmonic.h diff --git a/gsVisitorBiharmonicMixed.h b/src/gsVisitorBiharmonicMixed.h similarity index 100% rename from gsVisitorBiharmonicMixed.h rename to src/gsVisitorBiharmonicMixed.h diff --git a/gsVisitorElPoisson.h b/src/gsVisitorElPoisson.h similarity index 100% rename from gsVisitorElPoisson.h rename to src/gsVisitorElPoisson.h diff --git a/gsVisitorElUtils.h b/src/gsVisitorElUtils.h similarity index 60% rename from gsVisitorElUtils.h rename to src/gsVisitorElUtils.h index 336f159..a37b8a1 100644 --- a/gsVisitorElUtils.h +++ b/src/gsVisitorElUtils.h @@ -55,6 +55,36 @@ inline void symmetricIdentityTensor(gsMatrix & C, const gsMatrix & R) R(voigt(dim,i,0),voigt(dim,j,1))*R(voigt(dim,i,1),voigt(dim,j,0))); } +// Deviatoric matrix operator in 2D and 3D +template //template it with dimension +inline void deviatoricTensor(gsMatrix & C, const gsMatrix & R) +{ + short_t dim = R.cols(); + short_t dimTensor = dim*(dim+1)/2; + C.setZero(dimTensor,dimTensor); + + C.block(0,0,dim,dim).setConstant(-1.0/3.0); + C.block(0,0,dim,dim).diagonal().setConstant(2.0/3.0); + C.block(dim,dim,dimTensor-dim,dimTensor-dim).diagonal().setConstant(0.5); + + // if (dim == 2) + // { + // C <<2.0/3.0, -1.0/3.0, 0., + // -1.0/3.0, 2.0/3.0, 0., + // 0., 0., 1.0/2.0; + // } + // else if (dim == 3) + // { + // C << 2.0/3.0, -1.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0, + // -1.0/3.0, 2.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0, + // -1.0/3.0, -1.0/3.0, 2.0/3.0, 0.0, 0.0, 0.0, + // 0.0, 0.0, 0.0, 1.0/2.0, 0.0, 0.0, + // 0.0, 0.0, 0.0, 0.0, 1.0/2.0, 0.0, + // 0.0, 0.0, 0.0, 0.0, 0.0, 1.0/2.0; + // } +} + + // construct a fourth order matrix-trace tensor C based on two second order symmetric tensors R and S // C_ijkl = R_ij*S_kl in Voigt notation template @@ -80,6 +110,51 @@ inline void voigtStress(gsVector & Svec, const gsMatrix & S) Svec(i) = S(voigt(dim,i,0),voigt(dim,i,1)); } +// transform stress Svec in Voigt notation to stress tensor S +template +inline void tensorStress(const index_t dim, const gsMatrix & Svec, gsMatrix & S) +{ + if (dim == 2) + { + S.setZero(dim,dim); + S(0,0) = Svec(0,0); + S(1,1) = Svec(1,0); + S(0,1) = Svec(2,0); + S(1,0) = Svec(2,0); + } + else if (dim == 3) + { + S.setZero(dim,dim); + S(0,0) = Svec(0,0); + S(1,1) = Svec(1,0); + S(2,2) = Svec(2,0); + + S(0,1) = Svec(3,0); // tao_xy + S(1,0) = Svec(3,0); // tao_xy + + S(0,2) = Svec(4,0); // tao_yz + S(2,0) = Svec(4,0); // tao_yz + + S(1,2) = Svec(5,0); // tao_xz + S(2,1) = Svec(5,0); // tao_xz + } +} + +// transform strain tensor E to a vector in Voigt notation +template +inline void voigtStrain(gsVector & Evec, const gsMatrix & E) +{ + short_t dim = E.cols(); + //gsDebugVar(dim); + short_t dimTensor = dim*(dim+1)/2; + Evec.resize(dimTensor); + for (short i = 0; i < dimTensor; ++i) + if (voigt(dim,i,0) != voigt(dim,i,1)) + Evec(i) = E(voigt(dim,i,0),voigt(dim,i,1)) + E(voigt(dim,i,1),voigt(dim,i,0)); // off-diagonal terms + else + Evec(i) = E(voigt(dim,i,0),voigt(dim,i,1)); // diagonal terms +} + // auxiliary matrix B such that E:S = B*Svec in the weak form // (see Bernal, Calo, Collier, et. at., "PetIGA", ICCS 2013, p. 1610) template diff --git a/gsVisitorElasticityNeumann.h b/src/gsVisitorElasticityNeumann.h similarity index 100% rename from gsVisitorElasticityNeumann.h rename to src/gsVisitorElasticityNeumann.h diff --git a/gsVisitorLinearElasticity.h b/src/gsVisitorLinearElasticity.h similarity index 97% rename from gsVisitorLinearElasticity.h rename to src/gsVisitorLinearElasticity.h index 27940a3..df92441 100644 --- a/gsVisitorLinearElasticity.h +++ b/src/gsVisitorLinearElasticity.h @@ -30,8 +30,11 @@ class gsVisitorLinearElasticity public: gsVisitorLinearElasticity(const gsPde & pde_, gsSparseMatrix * elimMatrix = nullptr) - : dim(0), pde_ptr(static_cast*>(&pde_)), N_D(0), - elimMat(elimMatrix) + : + dim(0), + pde_ptr(static_cast*>(&pde_)), + N_D(0), + elimMat(elimMatrix) {} void initialize(const gsBasisRefs & basisRefs, @@ -103,7 +106,7 @@ class gsVisitorLinearElasticity { // stiffness matrix K = B_i^T * C * B_j; setB(B_i,I,physGrad.col(i)); - tempK = B_i.transpose() * C; + tempK = B_i.transpose() * C; //.reshapeCol(q,dim,dim) // loop over active basis functions (v_j) for (index_t j = 0; j < N_D; j++) { diff --git a/src/gsVisitorLinearElasticityMM.h b/src/gsVisitorLinearElasticityMM.h new file mode 100644 index 0000000..ef56e2b --- /dev/null +++ b/src/gsVisitorLinearElasticityMM.h @@ -0,0 +1,204 @@ +/** @file gsVisitorLinearElasticityMM.h + + @brief Visitor class for volumetric integration of the linear elasticity system. + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + O. Weeger (2012 - 2015, TU Kaiserslautern), + A.Shamanskiy (2016 - ...., TU Kaiserslautern) +*/ + +#pragma once + +#include +#include +#include +#include +#include + +#include +#include + +namespace gismo +{ + +template +class gsVisitorLinearElasticityMM //material matrix +{ +public: + + gsVisitorLinearElasticityMM(const gsPde & pde_,// const gsPieceWiseFunction &mmat, + const gsMaterialContainer & materials, + gsSparseMatrix * elimMatrix = nullptr) + : + pde_ptr(static_cast*>(&pde_)), + m_materials(materials), + elimMat(elimMatrix) + {} + + void initialize(const gsBasisRefs & basisRefs, + const index_t patchIndex, + const gsOptionList & options, + gsQuadRule & rule) + { + // parametric dimension of the first displacement component + dim = basisRefs.front().dim(); + // a quadrature rule is defined by the basis for the first displacement component. + rule = gsQuadrature::get(basisRefs.front(), options); + forceScaling = options.getReal("ForceScaling"); + localStiffening = options.getReal("LocalStiff"); + // saving necessary info + //--------------------------------------- + // T E = options.getReal("YoungsModulus"); + // T pr = options.getReal("PoissonsRatio"); + // m_materialMat = new gsElasticityMaterial(E,pr,dim); + //--------------------------------------- + I = gsMatrix::Identity(dim,dim); + // resize containers for global indices + globalIndices.resize(dim); + blockNumbers.resize(dim); + } + + inline void evaluate(const gsBasisRefs & basisRefs, + const gsGeometry & geo, + const gsMatrix & quNodes) + { + // store quadrature points of the element for geometry evaluation + md.points = quNodes; + // NEED_VALUE to get points in the physical domain for evaluation of the RHS + // NEED_MEASURE to get the Jacobian determinant values for integration + // NEED_GRAD_TRANSFORM to get the Jacobian matrix to transform gradient from the parametric to physical domain + md.flags = NEED_VALUE | NEED_MEASURE | NEED_GRAD_TRANSFORM; + // Compute image of the quadrature points plus gradient, jacobian and other necessary data + geo.computeMap(md); + // find local indices of the displacement basis functions active on the element + basisRefs.front().active_into(quNodes.col(0),localIndicesDisp); + N_D = localIndicesDisp.rows(); + // Evaluate displacement basis functions and their derivatives on the element + basisRefs.front().evalAllDers_into(quNodes,1,basisValuesDisp); + // Evaluate right-hand side at the image of the quadrature points + pde_ptr->rhs()->eval_into(md.values[0],forceValues); + + // Use the precomputed geometry data to evaluate the material matrix + gsMapData mdDeformed = md; + + gsMaterialData data; + gsMaterialBase * material = m_materials.piece(geo.id()); + material->precompute(md,mdDeformed,data); + material->eval_matrix_into(data,matValues); + + } + + inline void assemble(gsDomainIteratorWrapper & element, + const gsVector & quWeights) + { + // initialize local matrix and rhs + localMat.setZero(dim*N_D,dim*N_D); + localRhs.setZero(dim*N_D,1); + // Loop over the quadrature nodes + for (index_t q = 0; q < quWeights.rows(); ++q) + { + // Multiply quadrature weight by the geometry measure + const T weightForce = quWeights[q] * md.measure(q); + const T weightBody = quWeights[q] * math::pow(md.measure(q),1-localStiffening); + // Compute physical gradients of basis functions at q as a dim x numActiveFunction matrix + transformGradients(md,q,basisValuesDisp[1],physGrad); + C = matValues.reshapeCol(q,math::sqrt(matValues.rows()),math::sqrt(matValues.rows())); + // loop over active basis functions (v_j) + for (index_t i = 0; i < N_D; i++) + { + // stiffness matrix K = B_i^T * C * B_j; + setB(B_i,I,physGrad.col(i)); + tempK = B_i.transpose() * C; + // loop over active basis functions (v_j) + for (index_t j = 0; j < N_D; j++) + { + setB(B_j,I,physGrad.col(j)); + K = tempK * B_j; + for (short_t di = 0; di < dim; ++di) + for (short_t dj = 0; dj < dim; ++dj) + localMat(di*N_D+i,dj*N_D+j) += weightBody * K(di,dj); + } + } + // rhs contribution + for (short_t d = 0; d < dim; ++d) + localRhs.middleRows(d*N_D,N_D).noalias() += weightForce * forceScaling * forceValues(d,q) * basisValuesDisp[0].col(q) ; + } + } + + inline void localToGlobal(const int patchIndex, + const std::vector > & eliminatedDofs, + gsSparseSystem & system) + { + // computes global indices for displacement components + for (short_t d = 0; d < dim; ++d) + { + system.mapColIndices(localIndicesDisp, patchIndex, globalIndices[d], d); + blockNumbers.at(d) = d; + } + // push to global system + system.pushToRhs(localRhs,globalIndices,blockNumbers); + system.pushToMatrix(localMat,globalIndices,eliminatedDofs,blockNumbers,blockNumbers); + // push to the elimination system + if (elimMat != nullptr) + { + index_t globalI,globalElimJ; + index_t elimSize = 0; + for (short_t dJ = 0; dJ < dim; ++dJ) + { + for (short_t dI = 0; dI < dim; ++dI) + for (index_t i = 0; i < N_D; ++i) + if (system.colMapper(dI).is_free_index(globalIndices[dI].at(i))) + { + system.mapToGlobalRowIndex(localIndicesDisp.at(i),patchIndex,globalI,dI); + for (index_t j = 0; j < N_D; ++j) + if (!system.colMapper(dJ).is_free_index(globalIndices[dJ].at(j))) + { + globalElimJ = system.colMapper(dJ).global_to_bindex(globalIndices[dJ].at(j)); + elimMat->coeffRef(globalI,elimSize+globalElimJ) += localMat(N_D*dI+i,N_D*dJ+j); + } + } + elimSize += eliminatedDofs[dJ].rows(); + } + } + } + +protected: + // problem info + short_t dim; + const gsBasePde * pde_ptr; + // Lame coefficients and force scaling factor + T lambda, mu, forceScaling, localStiffening; + // geometry mapping + gsMapData md; + // local components of the global linear system + gsMatrix localMat; + gsMatrix localRhs; + // local indices (at the current patch) of the displacement basis functions active at the current element + gsMatrix localIndicesDisp; + // number of displacement basis functions active at the current element + index_t N_D; + // values and derivatives of displacement basis functions at quadrature points at the current element + // values are stored as a N_D x numQuadPoints matrix; not sure about derivatives, must be smth like N_D*dim x numQuadPoints + std::vector > basisValuesDisp; + // RHS values at quadrature points at the current element; stored as a dim x numQuadPoints matrix + gsMatrix forceValues; + /// Material handler + const gsMaterialContainer & m_materials; + // elimination matrix to efficiently change Dirichlet degrees of freedom + gsSparseMatrix * elimMat; + // all temporary matrices defined here for efficiency + gsMatrix C, Ctemp,physGrad, B_i, tempK, B_j, K, I; + // containers for global indices + std::vector< gsMatrix > globalIndices; + gsVector blockNumbers; + + gsMatrix matValues; +}; + +} // namespace gismo diff --git a/gsVisitorMass.h b/src/gsVisitorMass.h similarity index 100% rename from gsVisitorMass.h rename to src/gsVisitorMass.h diff --git a/gsVisitorMassElasticity.h b/src/gsVisitorMassElasticity.h similarity index 100% rename from gsVisitorMassElasticity.h rename to src/gsVisitorMassElasticity.h diff --git a/gsVisitorMixedLinearElasticity.h b/src/gsVisitorMixedLinearElasticity.h similarity index 98% rename from gsVisitorMixedLinearElasticity.h rename to src/gsVisitorMixedLinearElasticity.h index d4351f5..6948990 100644 --- a/gsVisitorMixedLinearElasticity.h +++ b/src/gsVisitorMixedLinearElasticity.h @@ -29,7 +29,11 @@ class gsVisitorMixedLinearElasticity { public: gsVisitorMixedLinearElasticity(const gsPde & pde_) - : dim(0), pde_ptr(static_cast*>(&pde_)), N_D(0) {} + : + dim(0), + pde_ptr(static_cast*>(&pde_)), + N_D(0) + {} void initialize(const gsBasisRefs & basisRefs, const index_t patchIndex, diff --git a/gsVisitorMixedNonLinearElasticity.h b/src/gsVisitorMixedNonLinearElasticity.h similarity index 100% rename from gsVisitorMixedNonLinearElasticity.h rename to src/gsVisitorMixedNonLinearElasticity.h diff --git a/gsVisitorMuscle.h b/src/gsVisitorMuscle.h similarity index 100% rename from gsVisitorMuscle.h rename to src/gsVisitorMuscle.h diff --git a/gsVisitorNavierStokes.h b/src/gsVisitorNavierStokes.h similarity index 100% rename from gsVisitorNavierStokes.h rename to src/gsVisitorNavierStokes.h diff --git a/gsVisitorNonLinearElasticity.h b/src/gsVisitorNonLinearElasticity.h similarity index 100% rename from gsVisitorNonLinearElasticity.h rename to src/gsVisitorNonLinearElasticity.h diff --git a/src/gsVisitorNonLinearElasticityMM.h b/src/gsVisitorNonLinearElasticityMM.h new file mode 100644 index 0000000..176063b --- /dev/null +++ b/src/gsVisitorNonLinearElasticityMM.h @@ -0,0 +1,284 @@ +/** @file gsVisitorNonLinearElasticityMM.h + + @brief Element visitor for nonlinear elasticity for 2D plain strain and 3D continua. + + This file is part of the G+Smo library. + + This Source Code Form is subject to the terms of the Mozilla Public + License, v. 2.0. If a copy of the MPL was not distributed with this + file, You can obtain one at http://mozilla.org/MPL/2.0/. + + Author(s): + O. Weeger (2012 - 2015, TU Kaiserslautern), + A.Shamanskiy (2016 - ...., TU Kaiserslautern) +*/ + +#pragma once + +#include +#include +#include +#include + +#include +#include + +namespace gismo +{ + +// template +// class gsAsDeformed : public gsFunction +// { +// public: +// gsAsDeformed(const gsFunction & undeformed_, +// const gsFunction & deformed_) +// : +// undeformed(undeformed_), +// deformed(deformed_) +// { +// } + +// short_t targetDim() const override { return undeformed.targetDim(); } +// short_t domainDim() const override { return undeformed.domainDim(); } + +// void eval_into(const gsMatrix & u, gsMatrix & result) const override +// { +// undeformed.eval_into(u,result); +// result += deformed.eval(u); +// } + +// void deriv_into(const gsMatrix & u, gsMatrix & result) const override +// { +// undeformed.deriv_into(u,result); +// result += deformed.deriv(u); +// } + +// void deriv2_into(const gsMatrix & u, gsMatrix & result) const override +// { +// undeformed.deriv2_into(u,result); +// result += deformed.deriv2(u); +// } + +// void evalAllDers_into(const gsMatrix & u, const int n, std::vector > & result, bool sameElement) const +// { +// std::vector> defResult; +// undeformed.evalAllDers_into(u,n,result,sameElement); +// deformed.evalAllDers_into(u,n,defResult,sameElement); +// for (index_t i = 0; i <= n; ++i) +// result[i] += defResult[i]; +// } + +// std::ostream &print(std::ostream &os) const +// { +// gsInfo<<"As deformed function\n"; +// undeformed.print(os); +// deformed.print(os); +// return os; +// } + +// GISMO_CLONE_FUNCTION(gsAsDeformed) + +// protected: +// const gsFunction & undeformed; +// const gsFunction & deformed; +// }; + + + +template +class gsVisitorNonLinearElasticityMM +{ +public: + gsVisitorNonLinearElasticityMM( const gsPde & pde_, + const gsMaterialContainer & materials, + const gsMultiPatch & displacement_) + : + pde_ptr(static_cast*>(&pde_)), + m_materials(materials), + displacement(displacement_) + { + } + + void initialize(const gsBasisRefs & basisRefs, + const index_t patchIndex, + const gsOptionList & options, + gsQuadRule & rule) + { + // parametric dimension of the first displacement component + dim = basisRefs.front().dim(); + // a quadrature rule is defined by the basis for the first displacement component. + rule = gsQuadrature::get(basisRefs.front(), options); + // saving necessary info + patch = patchIndex; + forceScaling = options.getReal("ForceScaling"); + localStiffening = options.getReal("LocalStiff"); + I = gsMatrix::Identity(dim,dim); + // resize containers for global indices + globalIndices.resize(dim); + blockNumbers.resize(dim); + } + + inline void evaluate(const gsBasisRefs & basisRefs, + const gsGeometry & geo, + const gsMatrix & quNodes) + { + GISMO_ASSERT((index_t)geo.id()==patch,"Geometry id ("<rhs()->eval_into(md.values[0],forceValues); + // store quadrature points of the element for displacement evaluation + mdDisplacement.points = quNodes; + // NEED_DERIV to compute deformation gradient + mdDisplacement.flags = NEED_DERIV; + // evaluate displacement gradient + displacement.patch(patch).computeMap(mdDisplacement); + + mdDeformed = mdDisplacement; + mdDeformed.values[1] += md.values[1]; + + // gsAsDeformed def(geo,displacement.patch(patch)); + + // // Construct a material evaluator for matrix with id geo.id() + // gsMaterialEvalSingle Seval(0,m_materials.piece(geo.id()),geo,def); + // gsMaterialEvalSingle Ceval(0,m_materials.piece(geo.id()),geo,def); + + // // The evaluator is single-piece, hence we use piece(0) + // Seval.eval_into(quNodes,vecValues); + // Ceval.eval_into(quNodes,matValues); + + // The lines below are faster than calling gsMaterialEval(Single), + // Since we re-use the geometric data and compute the parameter data only once + gsMaterialData data; + gsMaterialBase * material = m_materials.piece(geo.id()); + material->precompute(md,mdDeformed,data); + material->eval_deformation_gradient_into(data,defGradValues); // re-use the pre-computed ones + material->eval_stress_into(data,vecValues); + material->eval_matrix_into(data,matValues); + } + + inline void assemble(gsDomainIteratorWrapper & element, + const gsVector & quWeights) + { + // initialize local matrix and rhs + localMat.setZero(dim*N_D,dim*N_D); + localRhs.setZero(dim*N_D,1); + // loop over quadrature nodes + for (index_t q = 0; q < quWeights.rows(); ++q) + { + const T weightForce = quWeights[q] * md.measure(q); + // Compute physical gradients of basis functions at q as a dim x numActiveFunction matrix + transformGradients(md,q,basisValuesDisp[1],physGrad); + // deformation gradient + F = defGradValues.reshapeCol(q,dim,dim); + const T weightBody = quWeights[q] * pow(md.measure(q),-1.*localStiffening) * md.measure(q); + // Elasticity tensor + C = matValues.reshapeCol(q,math::sqrt(matValues.rows()),math::sqrt(matValues.rows())); + // Second Piola-Kirchhoff stress tensor + S = vecValues.reshapeCol(q,math::sqrt(vecValues.rows()),math::sqrt(vecValues.rows())); + for (index_t i = 0; i < N_D; i++) + { + // Material tangent K_tg_mat = B_i^T * C * B_j; + setB(B_i,F,physGrad.col(i)); + materialTangentTemp = B_i.transpose() * C; + // Geometric tangent K_tg_geo = gradB_i^T * S * gradB_j; + geometricTangentTemp = S * physGrad.col(i); + // loop over active basis functions (v_j) + for (index_t j = 0; j < N_D; j++) + { + setB(B_j,F,physGrad.col(j)); + + materialTangent = materialTangentTemp * B_j; + T geometricTangent = geometricTangentTemp.transpose() * physGrad.col(j); + // K_tg = K_tg_mat + I*K_tg_geo; + for (short_t d = 0; d < dim; ++d) + materialTangent(d,d) += geometricTangent; + + for (short_t di = 0; di < dim; ++di) + for (short_t dj = 0; dj < dim; ++dj) + localMat(di*N_D+i, dj*N_D+j) += weightBody * materialTangent(di,dj); + } + // Second Piola-Kirchhoff stress tensor as vector + voigtStress(Svec,S); + // rhs = -r = force - B*Svec, + localResidual = B_i.transpose() * Svec; + for (short_t d = 0; d < dim; d++) + localRhs(d*N_D+i) -= weightBody * localResidual(d); + } + // contribution of volumetric load function to residual/rhs + for (short_t d = 0; d < dim; ++d) + localRhs.middleRows(d*N_D,N_D).noalias() += weightForce * forceScaling * forceValues(d,q) * basisValuesDisp[0].col(q); + } + } + + inline void localToGlobal(const int patchIndex, + const std::vector > & eliminatedDofs, + gsSparseSystem & system) + { + // computes global indices for displacement components + for (short_t d = 0; d < dim; ++d) + { + system.mapColIndices(localIndicesDisp, patchIndex, globalIndices[d], d); + blockNumbers.at(d) = d; + } + // push to global system + system.pushToRhs(localRhs,globalIndices,blockNumbers); + system.pushToMatrix(localMat,globalIndices,eliminatedDofs,blockNumbers,blockNumbers); + } + +protected: + // problem info + short_t dim; + index_t patch; // current patch + const gsBasePde * pde_ptr; + index_t materialLaw; // (0: St. Venant-Kirchhoff, 1: ln neo-Hooke, 2: quad neo-Hooke) + // Lame coefficients and force scaling factor + T lambda, mu, forceScaling; + // geometry mapping + gsMapData md; + // local components of the global linear system + gsMatrix localMat; + gsMatrix localRhs; + // local indices (at the current patch) of the displacement basis functions active at the current element + gsMatrix localIndicesDisp; + // number of displacement basis functions active at the current element + index_t N_D; + // values and derivatives of displacement basis functions at quadrature points at the current element + // values are stored as a N_D x numQuadPoints matrix; not sure about derivatives, must be smth like N_D*dim x numQuadPoints + std::vector > basisValuesDisp; + // RHS values at quadrature points at the current element; stored as a dim x numQuadPoints matrix + gsMatrix forceValues; + /// Material handler + const gsMaterialContainer & m_materials; + // current displacement field + const gsMultiPatch & displacement; + // evaluation data of the current displacement field + gsMapData mdDisplacement; + gsMapData mdDeformed; + + // all temporary matrices defined here for efficiency + gsMatrix C, Ctemp, physGrad, physDispJac, F, RCG, E, S, RCGinv, B_i, materialTangentTemp, B_j, materialTangent, I; + gsVector geometricTangentTemp, Svec, localResidual; + T localStiffening; + // containers for global indices + std::vector< gsMatrix > globalIndices; + gsVector blockNumbers; + + gsMatrix matValues; + gsMatrix vecValues; + gsMatrix defGradValues; + +}; + +} // namespace gismo diff --git a/gsVisitorStokes.h b/src/gsVisitorStokes.h similarity index 100% rename from gsVisitorStokes.h rename to src/gsVisitorStokes.h diff --git a/gsVisitorThermo.h b/src/gsVisitorThermo.h similarity index 100% rename from gsVisitorThermo.h rename to src/gsVisitorThermo.h diff --git a/gsVisitorThermoBoundary.h b/src/gsVisitorThermoBoundary.h similarity index 100% rename from gsVisitorThermoBoundary.h rename to src/gsVisitorThermoBoundary.h diff --git a/gsWriteParaviewMultiPhysics.h b/src/gsWriteParaviewMultiPhysics.h similarity index 100% rename from gsWriteParaviewMultiPhysics.h rename to src/gsWriteParaviewMultiPhysics.h diff --git a/gsWriteParaviewMultiPhysics.hpp b/src/gsWriteParaviewMultiPhysics.hpp similarity index 85% rename from gsWriteParaviewMultiPhysics.hpp rename to src/gsWriteParaviewMultiPhysics.hpp index 5f0fb6b..ae54505 100644 --- a/gsWriteParaviewMultiPhysics.hpp +++ b/src/gsWriteParaviewMultiPhysics.hpp @@ -27,51 +27,6 @@ namespace gismo { - -//---------- START REPEATED from gsWriteParaview.hpp - -template -void writeSingleControlNet(const gsGeometry & Geo, - std::string const & fn) -{ - const short_t d = Geo.parDim(); - gsMesh msh; - Geo.controlNet(msh); - const short_t n = Geo.geoDim(); - if ( n == 1 ) - { - gsMatrix anch = Geo.basis().anchors(); - // Lift vertices at anchor positions - for (std::size_t i = 0; i!= msh.numVertices(); ++i) - { - msh.vertex(i)[d] = msh.vertex(i)[0]; - msh.vertex(i).topRows(d) = anch.col(i); - } - } - else if (n>3) - { - gsDebug<<"Writing 4th coordinate\n"; - const gsMatrix & cp = Geo.coefs(); - gsWriteParaviewPoints(cp.transpose(), fn ); - return; - } - - gsWriteParaview(msh, fn, false); -} - -template -void writeSingleCompMesh(const gsBasis & basis, const gsGeometry & Geo, - std::string const & fn, unsigned resolution) -{ - gsMesh msh(basis, resolution); - Geo.evaluateMesh(msh); - gsWriteParaview(msh, fn, false); -} - - -//---------- END REPEATED from gsWriteParaview.hpp - - template void gsWriteParaviewMultiPhysics(std::map*> fields, std::string const & fn, diff --git a/gsWriteParaviewMultiPhysics_.cpp b/src/gsWriteParaviewMultiPhysics_.cpp similarity index 100% rename from gsWriteParaviewMultiPhysics_.cpp rename to src/gsWriteParaviewMultiPhysics_.cpp