From f4df77b29cf02d8631b5f4170826b4b0a13ed831 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Thu, 29 Feb 2024 17:00:44 +0100 Subject: [PATCH 1/4] temporary fix, using two `auto` functions registered to the expression assembler, and a switch to decide which one to use --- src/gsThinShellAssembler.h | 1 + src/gsThinShellAssembler.hpp | 67 +++++++++++++++++++++++++++--------- 2 files changed, 52 insertions(+), 16 deletions(-) diff --git a/src/gsThinShellAssembler.h b/src/gsThinShellAssembler.h index ec162d7..2724bc5 100644 --- a/src/gsThinShellAssembler.h +++ b/src/gsThinShellAssembler.h @@ -671,6 +671,7 @@ class gsThinShellAssembler : public gsThinShellAssemblerBase gsSparseMatrix m_mass; const gsFunction * m_forceFun; + bool m_parametricForce; const gsFunction * m_foundFun; const gsFunction * m_pressFun; diff --git a/src/gsThinShellAssembler.hpp b/src/gsThinShellAssembler.hpp index d691555..80a2b61 100644 --- a/src/gsThinShellAssembler.hpp +++ b/src/gsThinShellAssembler.hpp @@ -50,6 +50,9 @@ gsThinShellAssembler::gsThinShellAssembler(const gsMultiPatch m_forceFun(&surface_force), m_materialMatrices(materialMatrices) { + // surface forces defined in the parametric domain (ONLY WORKS FOR 3D) + m_parametricForce = (surface_force.domainDim()==2 && d==3); + this->_defaultOptions(); this->_getOptions(); this->_initialize(); @@ -75,6 +78,9 @@ gsThinShellAssembler::gsThinShellAssembler(const gsMultiPatch for (size_t p=0; p!=m_patches.nPatches(); p++) m_materialMatrices.set(p,materialMatrix); + // surface forces defined in the parametric domain (ONLY WORKS FOR 3D) + m_parametricForce = (surface_force.domainDim()==2 && d==3); + this->_defaultOptions(); this->_getOptions(); this->_initialize(); @@ -1442,6 +1448,8 @@ ThinShellAssemblerStatus gsThinShellAssembler::assembleMass(const else m_assembler.assemble(mm0.val()*(m_space.rowSum())*meas(m_ori)); + m_mass = m_assembler.matrix(); + /* // assemble system if (!lumped) { @@ -1545,8 +1553,8 @@ gsThinShellAssembler::assemble_impl() space m_space = m_assembler.trialSpace(0); - auto m_force = m_assembler.getCoeff(*m_forceFun, m_ori); - + auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain + auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ); auto m_Ef_der = -( deriv2(m_space,sn(m_def).normalized().tr() ) + deriv2(m_def,var1(m_space,m_def) ) ) * reshape(m_m2,3,3); @@ -1573,10 +1581,13 @@ gsThinShellAssembler::assemble_impl() + m_M_der * m_Ef_der.tr() ) * meas(m_ori) - , - m_space * m_force * meas(m_ori) ); + if (m_parametricForce) // Assemble the force defined in the parameter domain + m_assembler.assemble(m_space * m_parforce * meas(m_ori)); + else // Assemble the force defined in the physical domain + m_assembler.assemble(m_space * m_physforce * meas(m_ori)); + this->_assembleWeakBCs(); this->_assembleWeakBCs(); this->_assembleWeakIfc(); @@ -1622,7 +1633,8 @@ gsThinShellAssembler::assemble_impl() auto mmA = m_assembler.getCoeff(m_mmA); space m_space = m_assembler.trialSpace(0); - auto m_force = m_assembler.getCoeff(*m_forceFun, m_ori); + auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain + auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain auto jacG = jac(m_def); auto m_Em_der = flat( jacG.tr() * jac(m_space) ) ; //[checked] @@ -1645,10 +1657,13 @@ gsThinShellAssembler::assemble_impl() ( m_N_der * m_Em_der.tr() ) * meas(m_ori) - , - m_space * m_force * meas(m_ori) ); + if (m_parametricForce) // Assemble the force defined in the parameter domain + m_assembler.assemble(m_space * m_parforce * meas(m_ori)); + else // Assemble the force defined in the physical domain + m_assembler.assemble(m_space * m_physforce * meas(m_ori)); + this->_assembleWeakBCs(); this->_assembleWeakBCs(); this->_assembleWeakIfc(); @@ -1983,7 +1998,8 @@ gsThinShellAssembler::assembleVector_impl(const gsFunctionSet auto m_m2 = m_assembler.getCoeff(mult2t); space m_space = m_assembler.trialSpace(0); - auto m_force = m_assembler.getCoeff(*m_forceFun, m_ori); + auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain + auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain if (homogenize) this->homogenizeDirichlet(); @@ -2000,10 +2016,19 @@ gsThinShellAssembler::assembleVector_impl(const gsFunctionSet if (m_foundInd) this->_assembleFoundation(*m_foundFun,deformed); if (m_pressInd) this->_assemblePressure(*m_pressFun,deformed); - // Assemble vector - m_assembler.assemble(m_space * m_force * meas(m_ori) - - ( ( m_N * m_Em_der.tr() + m_M * m_Ef_der.tr() ) * meas(m_ori) ).tr() - ); + // Assemble vector + ////// External force + if (m_parametricForce) // Assemble the force defined in the parameter domain + m_assembler.assemble(m_space * m_parforce * meas(m_ori)); + else // Assemble the force defined in the physical domain + m_assembler.assemble(m_space * m_physforce * meas(m_ori)); + + ////// Internal force + m_assembler.assemble( + ( + - ( ( m_N * m_Em_der.tr() + m_M * m_Ef_der.tr() ) * meas(m_ori) ).tr() + ) + ); this->_assembleWeakBCs(deformed); this->_assembleWeakIfc(deformed); @@ -2045,7 +2070,8 @@ gsThinShellAssembler::assembleVector_impl(const gsFunctionSet auto S0 = m_assembler.getCoeff(m_S0); space m_space = m_assembler.trialSpace(0); - auto m_force = m_assembler.getCoeff(*m_forceFun, m_ori); + auto m_physforce = m_assembler.getCoeff(*m_forceFun,m_ori); // force defined in physical domain + auto m_parforce = m_assembler.getCoeff(*m_forceFun); // force defined in parametric domain if (homogenize) this->homogenizeDirichlet(); else this->_assembleDirichlet(); @@ -2059,9 +2085,18 @@ gsThinShellAssembler::assembleVector_impl(const gsFunctionSet if (m_pressInd) this->_assemblePressure(*m_pressFun,deformed); // Assemble vector - m_assembler.assemble(m_space * m_force * meas(m_ori) - - ( ( m_N * m_Em_der.tr() ) * meas(m_ori) ).tr() - ); + ////// External force + if (m_parametricForce) // Assemble the force defined in the parameter domain + m_assembler.assemble(m_space * m_parforce * meas(m_ori)); + else // Assemble the force defined in the physical domain + m_assembler.assemble(m_space * m_physforce * meas(m_ori)); + + ////// Internal force + m_assembler.assemble( + ( + - ( ( m_N * m_Em_der.tr() ) * meas(m_ori) ).tr() + ) + ); this->_assembleWeakBCs(deformed); this->_assembleWeakIfc(deformed); From 57275857fd4608ca20afa28c80fcb07f1a5c5bab Mon Sep 17 00:00:00 2001 From: Crazy-Rich-Meghan Date: Thu, 11 Apr 2024 11:42:29 +0200 Subject: [PATCH 2/4] Constructor for the ThinShellAssembler which takes gsFunctionSet as a parameter --- src/gsMaterialMatrixLinear.h | 2 +- src/gsMaterialMatrixNonlinear.h | 2 +- src/gsThinShellAssembler.h | 6 +++--- src/gsThinShellAssembler.hpp | 4 ++-- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/gsMaterialMatrixLinear.h b/src/gsMaterialMatrixLinear.h index 3fed637..179b5b6 100644 --- a/src/gsMaterialMatrixLinear.h +++ b/src/gsMaterialMatrixLinear.h @@ -42,7 +42,7 @@ class gsMaterialMatrixLinear : public gsMaterialMatrixBaseDim typedef T Scalar_t; - GISMO_OVERRIDE_CLONE_FUNCTION(gsMaterialMatrixLinear) + GISMO_OVERRIDE_CLONE_FUNCTION(gsMaterialMatrixLinear); using Base = gsMaterialMatrixBaseDim; diff --git a/src/gsMaterialMatrixNonlinear.h b/src/gsMaterialMatrixNonlinear.h index f8e0831..d192ea7 100644 --- a/src/gsMaterialMatrixNonlinear.h +++ b/src/gsMaterialMatrixNonlinear.h @@ -48,7 +48,7 @@ class gsMaterialMatrixNonlinear : public gsMaterialMatrixBaseDim { public: - GISMO_OVERRIDE_CLONE_FUNCTION(gsMaterialMatrixNonlinear) + GISMO_OVERRIDE_CLONE_FUNCTION(gsMaterialMatrixNonlinear); using Base = gsMaterialMatrixBaseDim; diff --git a/src/gsThinShellAssembler.h b/src/gsThinShellAssembler.h index 2724bc5..52b92a3 100644 --- a/src/gsThinShellAssembler.h +++ b/src/gsThinShellAssembler.h @@ -92,7 +92,7 @@ class gsThinShellAssembler : public gsThinShellAssemblerBase gsThinShellAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, const gsBoundaryConditions & bconditions, - const gsFunction & surface_force, + const gsFunctionSet & surface_force, const gsMaterialMatrixContainer & materialmatrices); /** @@ -107,7 +107,7 @@ class gsThinShellAssembler : public gsThinShellAssemblerBase gsThinShellAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, const gsBoundaryConditions & bconditions, - const gsFunction & surface_force, + const gsFunctionSet & surface_force, gsMaterialMatrixBase * materialmatrix); /** @@ -670,7 +670,7 @@ class gsThinShellAssembler : public gsThinShellAssemblerBase gsSparseMatrix m_mass; - const gsFunction * m_forceFun; + const gsFunctionSet * m_forceFun; bool m_parametricForce; const gsFunction * m_foundFun; const gsFunction * m_pressFun; diff --git a/src/gsThinShellAssembler.hpp b/src/gsThinShellAssembler.hpp index 80a2b61..2a4ae39 100644 --- a/src/gsThinShellAssembler.hpp +++ b/src/gsThinShellAssembler.hpp @@ -39,7 +39,7 @@ template gsThinShellAssembler::gsThinShellAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, const gsBoundaryConditions & bconditions, - const gsFunction & surface_force, + const gsFunctionSet & surface_force, const gsMaterialMatrixContainer & materialMatrices ) : @@ -62,7 +62,7 @@ template gsThinShellAssembler::gsThinShellAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, const gsBoundaryConditions & bconditions, - const gsFunction & surface_force, + const gsFunctionSet & surface_force, gsMaterialMatrixBase * materialMatrix ) : From 6a3c519a3b6f7ea48846062202ddc77d3f0566f0 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Wed, 24 Apr 2024 10:45:03 +0200 Subject: [PATCH 3/4] add references to foundation implementations --- src/gsThinShellAssembler.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/gsThinShellAssembler.hpp b/src/gsThinShellAssembler.hpp index 2a4ae39..ada469e 100644 --- a/src/gsThinShellAssembler.hpp +++ b/src/gsThinShellAssembler.hpp @@ -407,6 +407,7 @@ void gsThinShellAssembler::_assemblePressure(const gsFunction _assemblePressure_impl(pressFun,deformed); } +// assembles eq 3.26 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78 template template typename std::enable_if<(_d==3) && matrix, void>::type @@ -426,6 +427,7 @@ gsThinShellAssembler::_assemblePressure_impl(const gsFunction ); } +// assembles eq 3.25 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78 template template typename std::enable_if<(_d==3) && !_matrix, void>::type @@ -469,6 +471,7 @@ gsThinShellAssembler::_assembleFoundation_impl(const gsFunction template typename std::enable_if<(_d==3) && !_matrix, void>::type @@ -502,6 +505,7 @@ void gsThinShellAssembler::_assembleFoundation(const gsFunction(foundFun,deformed); } +// assembles eq 3.28 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78 template template typename std::enable_if<(_d==3) && matrix, void>::type From 0b78123c5120f37548249f9f0e9c7ee3baaf9c97 Mon Sep 17 00:00:00 2001 From: Hugo Verhelst Date: Tue, 7 May 2024 15:59:22 +0200 Subject: [PATCH 4/4] foundation bug --- src/gsThinShellAssembler.hpp | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/src/gsThinShellAssembler.hpp b/src/gsThinShellAssembler.hpp index ada469e..64db781 100644 --- a/src/gsThinShellAssembler.hpp +++ b/src/gsThinShellAssembler.hpp @@ -467,15 +467,6 @@ template template typename std::enable_if<(_d==3) && _matrix, void>::type gsThinShellAssembler::_assembleFoundation_impl(const gsFunction & foundFun) -{ - // No matrix contribution for the linear case -} - -// assembles eq 3.27 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78 -template -template -typename std::enable_if<(_d==3) && !_matrix, void>::type -gsThinShellAssembler::_assembleFoundation_impl(const gsFunction & foundFun) { geometryMap m_ori = m_assembler.getMap(m_patches); @@ -489,6 +480,15 @@ gsThinShellAssembler::_assembleFoundation_impl(const gsFunction +template +typename std::enable_if<(_d==3) && !_matrix, void>::type +gsThinShellAssembler::_assembleFoundation_impl(const gsFunction & foundFun) +{ + // No rhs contribution for the linear case +} + template template typename std::enable_if::type @@ -1733,10 +1733,6 @@ gsThinShellAssembler::assembleMatrix_impl(const gsFunctionSet this->homogenizeDirichlet(); - gsVector pt(2); - pt.setConstant(0.25); - gsExprEvaluator ev(m_assembler); - auto m_N = S0.tr(); auto m_Em_der = flat( jac(m_def).tr() * jac(m_space) ) ; //[checked] auto m_Em_der2 = flatdot( jac(m_space),jac(m_space).tr(), m_N ); //[checked]