Add support for material classes

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} $<TARGET_OBJECTS:${PROJECT_NAME}>
   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" )

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 <gismo.h>
+#include <gsElasticity/gsElasticityAssembler.h>
+#include <gsElasticity/gsIterative.h>
+#include <gsElasticity/gsWriteParaviewMultiPhysics.h>
+#include <gsElasticity/gsMaterialBase.h>
+#include <gsElasticity/gsLinearMaterial.h>
+
+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<real_t> materialMat(youngsModulus,poissonsRatio,3);
+
+    // creating assembler
+    gsElasticityAssembler<real_t> assembler(geometry,basis,bcInfo,f,&materialMat);
+    // gsElasticityAssembler<real_t> 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<real_t> 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<std::string,const gsField<> *> 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;
+}

examples/cooks_nonLinElast_2D.cpp

@@ -7,6 +7,9 @@
 #include <gsElasticity/gsElasticityAssembler.h>
 #include <gsElasticity/gsIterative.h>
 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
+#include <gsElasticity/gsMaterialBase.h>
+#include <gsElasticity/gsLinearMaterial.h>
+#include <gsElasticity/gsNeoHookeLogMaterial.h>
 
 using namespace gismo;
 
@@ -67,11 +70,16 @@ int main(int argc, char* argv[]){
                   // Solving //
     //=============================================//
 
+    gsNeoHookeLogMaterial<real_t> materialMat(youngsModulus,poissonsRatio,2);
+    // gsLinearMaterial<real_t> materialMat(youngsModulus,poissonsRatio);
+
     // creating assembler
-    gsElasticityAssembler<real_t> assembler(geometry,basisDisplacement,bcInfo,g);
+    gsElasticityAssembler<real_t> assembler(geometry,basisDisplacement,bcInfo,g,&materialMat);
+    // gsElasticityAssembler<real_t> 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

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 <gismo.h>
+#include <gsElasticity/gsMaterialEval.h>
+#include <gsElasticity/gsLinearMaterial.h>
+#include <gsUtils/gsStopwatch.h>
+
+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<real_t> SvK(E_modulus, PoissonRatio, 3);
+    gsStopwatch time;
+    // ============================================================
+    time.restart();
+    // Using precomputation of material matrices
+    // Precompute the material matrix
+    gsMaterialData<real_t> 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 "<<time.stop()<<" s\n";
+
+    gsInfo<<"Deformation gradient\n";
+    gsInfo<<Fres<<"\n";
+
+    gsInfo<<"Strain\n";
+    gsInfo<<Eres<<"\n";
+
+    gsInfo<<"Stress\n";
+    gsInfo<<Sres<<"\n";
+
+    gsInfo<<"Material matrix\n";
+    gsInfo<<Cres<<"\n";
+
+    // ============================================================
+    // Using gsMaterialEval
+    time.restart();
+    gsMaterialEval<real_t, gsMaterialOutput::F, true> SvK_F(&SvK, mp, mp_def);
+    gsMaterialEval<real_t, gsMaterialOutput::E, true> SvK_E(&SvK, mp, mp_def);
+    gsMaterialEval<real_t, gsMaterialOutput::S, true> SvK_S(&SvK, mp, mp_def);
+    gsMaterialEval<real_t, gsMaterialOutput::C, true> 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 "<<time.stop()<<" s\n";
+
+    // gsInfo<<"Deformation gradient\n";
+    // gsInfo<<Fres<<"\n";
+
+    gsInfo<<"Strain\n";
+    gsInfo<<Eres<<"\n";
+
+    gsInfo<<"Stress\n";
+    gsInfo<<Sres<<"\n";
+
+    gsInfo<<"Material matrix\n";
+    gsInfo<<Cres<<"\n";
+
+    return 0;
+}

examples/plateWithHole_linElast_2D.cpp

@@ -6,6 +6,8 @@
 #include <gismo.h>
 #include <gsElasticity/gsElasticityAssembler.h>
 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
+#include <gsElasticity/gsMaterialBase.h>
+#include <gsElasticity/gsLinearMaterial.h>
 
 using namespace gismo;
 
@@ -75,8 +77,11 @@ int main(int argc, char* argv[]){
               // Assembling & solving //
     //=============================================//
 
+    gsLinearMaterial<real_t> materialMat(youngsModulus,poissonsRatio,2);
+
     // creating assembler
-    gsElasticityAssembler<real_t> assembler(geometry,basis,bcInfo,g);
+    // gsElasticityAssembler<real_t> assembler(geometry,basis,bcInfo,g);//,materialMat);
+    gsElasticityAssembler<real_t> assembler(geometry,basis,bcInfo,g,&materialMat);
     assembler.options().setReal("YoungsModulus",youngsModulus);
     assembler.options().setReal("PoissonsRatio",poissonsRatio);
     gsInfo<<"Assembling...\n";

examples/terrific_linElast_3D.cpp

@@ -6,7 +6,8 @@
 #include <gismo.h>
 #include <gsElasticity/gsElasticityAssembler.h>
 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
-#include <gsElasticity/gsGeoUtils.h>
+#include <gsElasticity/gsMaterialBase.h>
+#include <gsElasticity/gsLinearMaterial.h>
 
 using namespace gismo;
 
@@ -72,8 +73,12 @@ int main(int argc, char* argv[])
               // Assembling & solving //
     //=============================================//
 
+    gsLinearMaterial<real_t> materialMat(youngsModulus,poissonsRatio,3);
+
     // creating assembler
-    gsElasticityAssembler<real_t> assembler(geometry,basis,bcInfo,f);
+    // gsElasticityAssembler<real_t> assembler(geometry,basis,bcInfo,g);//,materialMat);
+    gsElasticityAssembler<real_t> assembler(geometry,basis,bcInfo,f,&materialMat);
+
     assembler.options().setReal("YoungsModulus",youngsModulus);
     assembler.options().setReal("PoissonsRatio",poissonsRatio);
     assembler.options().setInt("DirichletValues",dirichlet::l2Projection);

examples/terrific_nonLinElast_3D.cpp

@@ -6,6 +6,7 @@
 #include <gismo.h>
 #include <gsElasticity/gsElasticityAssembler.h>
 #include <gsElasticity/gsIterative.h>
+#include <gsElasticity/gsLinearMaterial.h>
 #include <gsElasticity/gsWriteParaviewMultiPhysics.h>
 
 using namespace gismo;
@@ -74,6 +75,8 @@ int main(int argc, char* argv[]){
                   // Solving //
     //=============================================//
 
+    gsLinearMaterial<real_t> materialMat(youngsModulus,poissonsRatio,3);
+
     // creating assembler
     gsElasticityAssembler<real_t> assembler(geometry,basis,bcInfo,f);
     assembler.options().setReal("YoungsModulus",youngsModulus);

gsALE.hpp → src/gsALE.hpp

@@ -42,7 +42,7 @@ gsALE<T>::gsALE(gsMultiPatch<T> & geometry, const gsMultiPatch<T> & 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<T> rhs(gsVector<>::Zero(geometry.parDim()),geometry.parDim());
+    gsConstantFunction<T> rhs(gsVector<T>::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)

gsBaseAssembler.h → src/gsBaseAssembler.h

@@ -92,7 +92,7 @@ class gsBaseAssembler : public gsAssembler<T>
 
     virtual void setRHS(const gsMatrix<T> & rhs) {m_system.rhs() = rhs;}
 
-    virtual void setMatrix(const gsSparseMatrix<T> & matrix) {m_system.matrix() = matrix;}
+    virtual void setMatrix(const gsSparseMatrix<T> & mat) {m_system.matrix() = mat;}
 
 protected:
     using gsAssembler<T>::m_pde_ptr;