Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Flexibility for parametric dimension of loads (Closes Issue #20) #21

Draft
wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Flexibility for parametric dimension of loads (Closes Issue #20) #21

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f4df77b
temporary fix, using two `auto` functions registered to the expressio…
hverhelst Feb 29, 2024
5727585
Constructor for the ThinShellAssembler which takes gsFunctionSet as a…
Crazy-Rich-Meghan Apr 11, 2024
6a3c519
add references to foundation implementations
hverhelst Apr 24, 2024
0b78123
foundation bug
hverhelst May 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/gsMaterialMatrixLinear.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ class gsMaterialMatrixLinear : public gsMaterialMatrixBaseDim<dim,T>

typedef T Scalar_t;

GISMO_OVERRIDE_CLONE_FUNCTION(gsMaterialMatrixLinear)
GISMO_OVERRIDE_CLONE_FUNCTION(gsMaterialMatrixLinear);

using Base = gsMaterialMatrixBaseDim<dim,T>;

Expand Down
2 changes: 1 addition & 1 deletion src/gsMaterialMatrixNonlinear.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ class gsMaterialMatrixNonlinear : public gsMaterialMatrixBaseDim<dim,T>
{
public:

GISMO_OVERRIDE_CLONE_FUNCTION(gsMaterialMatrixNonlinear)
GISMO_OVERRIDE_CLONE_FUNCTION(gsMaterialMatrixNonlinear);

using Base = gsMaterialMatrixBaseDim<dim,T>;

Expand Down
7 changes: 4 additions & 3 deletions src/gsThinShellAssembler.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ class gsThinShellAssembler : public gsThinShellAssemblerBase<T>
gsThinShellAssembler(const gsMultiPatch<T> & patches,
const gsMultiBasis<T> & basis,
const gsBoundaryConditions<T> & bconditions,
const gsFunction<T> & surface_force,
const gsFunctionSet<T> & surface_force,
const gsMaterialMatrixContainer<T> & materialmatrices);

/**
Expand All @@ -107,7 +107,7 @@ class gsThinShellAssembler : public gsThinShellAssemblerBase<T>
gsThinShellAssembler(const gsMultiPatch<T> & patches,
const gsMultiBasis<T> & basis,
const gsBoundaryConditions<T> & bconditions,
const gsFunction<T> & surface_force,
const gsFunctionSet<T> & surface_force,
gsMaterialMatrixBase<T> * materialmatrix);

/**
Expand Down Expand Up @@ -670,7 +670,8 @@ class gsThinShellAssembler : public gsThinShellAssemblerBase<T>

gsSparseMatrix<T> m_mass;

const gsFunction<T> * m_forceFun;
const gsFunctionSet<T> * m_forceFun;
bool m_parametricForce;
const gsFunction<T> * m_foundFun;
const gsFunction<T> * m_pressFun;

Expand Down
95 changes: 65 additions & 30 deletions src/gsThinShellAssembler.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ template<short_t d, class T, bool bending>
gsThinShellAssembler<d, T, bending>::gsThinShellAssembler(const gsMultiPatch<T> & patches,
const gsMultiBasis<T> & basis,
const gsBoundaryConditions<T> & bconditions,
const gsFunction<T> & surface_force,
const gsFunctionSet<T> & surface_force,
const gsMaterialMatrixContainer<T> & materialMatrices
)
:
Expand All @@ -50,6 +50,9 @@ gsThinShellAssembler<d, T, bending>::gsThinShellAssembler(const gsMultiPatch<T>
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();
Expand All @@ -59,7 +62,7 @@ template<short_t d, class T, bool bending>
gsThinShellAssembler<d, T, bending>::gsThinShellAssembler(const gsMultiPatch<T> & patches,
const gsMultiBasis<T> & basis,
const gsBoundaryConditions<T> & bconditions,
const gsFunction<T> & surface_force,
const gsFunctionSet<T> & surface_force,
gsMaterialMatrixBase<T> * materialMatrix
)
:
Expand All @@ -75,6 +78,9 @@ gsThinShellAssembler<d, T, bending>::gsThinShellAssembler(const gsMultiPatch<T>
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();
Expand Down Expand Up @@ -401,6 +407,7 @@ void gsThinShellAssembler<d, T, bending>::_assemblePressure(const gsFunction<T>
_assemblePressure_impl<d,_matrix>(pressFun,deformed);
}

// assembles eq 3.26 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78
template <short_t d, class T, bool bending>
template <short_t _d, bool matrix>
typename std::enable_if<(_d==3) && matrix, void>::type
Expand All @@ -420,6 +427,7 @@ gsThinShellAssembler<d, T, bending>::_assemblePressure_impl(const gsFunction<T>
);
}

// assembles eq 3.25 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78
template <short_t d, class T, bool bending>
template <short_t _d, bool _matrix>
typename std::enable_if<(_d==3) && !_matrix, void>::type
Expand Down Expand Up @@ -459,14 +467,6 @@ template <short_t d, class T, bool bending>
template <short_t _d, bool _matrix>
typename std::enable_if<(_d==3) && _matrix, void>::type
gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & foundFun)
{
// No matrix contribution for the linear case
}

template <short_t d, class T, bool bending>
template <short_t _d, bool _matrix>
typename std::enable_if<(_d==3) && !_matrix, void>::type
gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & foundFun)
{
geometryMap m_ori = m_assembler.getMap(m_patches);

Expand All @@ -480,6 +480,15 @@ gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T
);
}

// assembles eq 3.27 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78
template <short_t d, class T, bool bending>
template <short_t _d, bool _matrix>
typename std::enable_if<(_d==3) && !_matrix, void>::type
gsThinShellAssembler<d, T, bending>::_assembleFoundation_impl(const gsFunction<T> & foundFun)
{
// No rhs contribution for the linear case
}

template <short_t d, class T, bool bending>
template <short_t _d, bool _matrix>
typename std::enable_if<!(_d==3), void>::type
Expand All @@ -496,6 +505,7 @@ void gsThinShellAssembler<d, T, bending>::_assembleFoundation(const gsFunction<T
_assembleFoundation_impl<d,_matrix>(foundFun,deformed);
}

// assembles eq 3.28 from http://resolver.tudelft.nl/uuid:56c0cc91-643d-4817-9702-93fedce5fd78
template <short_t d, class T, bool bending>
template <short_t _d, bool matrix>
typename std::enable_if<(_d==3) && matrix, void>::type
Expand Down Expand Up @@ -1442,6 +1452,8 @@ ThinShellAssemblerStatus gsThinShellAssembler<d, T, bending>::assembleMass(const
else
m_assembler.assemble(mm0.val()*(m_space.rowSum())*meas(m_ori));

m_mass = m_assembler.matrix();

/* // assemble system
if (!lumped)
{
Expand Down Expand Up @@ -1545,8 +1557,8 @@ gsThinShellAssembler<d, T, bending>::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);
Expand All @@ -1573,10 +1585,13 @@ gsThinShellAssembler<d, T, bending>::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<true>();
this->_assembleWeakBCs<false>();
this->_assembleWeakIfc<true>();
Expand Down Expand Up @@ -1622,7 +1637,8 @@ gsThinShellAssembler<d, T, bending>::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]
Expand All @@ -1645,10 +1661,13 @@ gsThinShellAssembler<d, T, bending>::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<true>();
this->_assembleWeakBCs<false>();
this->_assembleWeakIfc<true>();
Expand Down Expand Up @@ -1714,10 +1733,6 @@ gsThinShellAssembler<d, T, bending>::assembleMatrix_impl(const gsFunctionSet<T>

this->homogenizeDirichlet();

gsVector<T> pt(2);
pt.setConstant(0.25);
gsExprEvaluator<T> 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]
Expand Down Expand Up @@ -1983,7 +1998,8 @@ gsThinShellAssembler<d, T, bending>::assembleVector_impl(const gsFunctionSet<T>
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();
Expand All @@ -2000,10 +2016,19 @@ gsThinShellAssembler<d, T, bending>::assembleVector_impl(const gsFunctionSet<T>
if (m_foundInd) this->_assembleFoundation<false>(*m_foundFun,deformed);
if (m_pressInd) this->_assemblePressure<false>(*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<false>(deformed);
this->_assembleWeakIfc<false>(deformed);
Expand Down Expand Up @@ -2045,7 +2070,8 @@ gsThinShellAssembler<d, T, bending>::assembleVector_impl(const gsFunctionSet<T>
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();
Expand All @@ -2059,9 +2085,18 @@ gsThinShellAssembler<d, T, bending>::assembleVector_impl(const gsFunctionSet<T>
if (m_pressInd) this->_assemblePressure<false>(*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<false>(deformed);
this->_assembleWeakIfc<false>(deformed);
Expand Down