diff --git a/src/gsMaterialMatrixLinear.h b/src/gsMaterialMatrixLinear.h index 9095217..80b71c4 100644 --- a/src/gsMaterialMatrixLinear.h +++ b/src/gsMaterialMatrixLinear.h @@ -40,7 +40,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 5683f11..39d9a21 100644 --- a/src/gsMaterialMatrixNonlinear.h +++ b/src/gsMaterialMatrixNonlinear.h @@ -47,7 +47,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 415ffb0..61f2d37 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); gsThinShellAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, @@ -675,7 +675,8 @@ 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 c98ad23..b575ce2 100644 --- a/src/gsThinShellAssembler.hpp +++ b/src/gsThinShellAssembler.hpp @@ -34,7 +34,7 @@ template gsThinShellAssembler::gsThinShellAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, const gsBoundaryConditions & bconditions, - const gsFunction & surface_force, + const gsFunctionSet & surface_force, const gsMaterialMatrixContainer & materialMatrices ) : @@ -45,6 +45,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(); @@ -54,7 +57,7 @@ template gsThinShellAssembler::gsThinShellAssembler(const gsMultiPatch & patches, const gsMultiBasis & basis, const gsBoundaryConditions & bconditions, - const gsFunction & surface_force, + const gsFunctionSet & surface_force, typename gsMaterialMatrixBase::uPtr & materialMatrix ) : @@ -83,6 +86,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(); @@ -409,6 +415,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 @@ -428,6 +435,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 @@ -467,14 +475,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 -} - -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); @@ -488,6 +488,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 @@ -504,6 +513,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 @@ -1450,6 +1460,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) { @@ -1553,8 +1565,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); @@ -1581,10 +1593,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(); @@ -1630,7 +1645,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] @@ -1653,10 +1669,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(); @@ -1722,10 +1741,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] @@ -1991,7 +2006,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(); @@ -2008,10 +2024,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); @@ -2053,7 +2078,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(); @@ -2067,9 +2093,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);