Skip to content

Commit

Permalink
clean up
Browse files Browse the repository at this point in the history
  • Loading branch information
@mohd-afeef-badri
mohd-afeef-badri committed Mar 4, 2025
1 parent e7ca372 commit 51c9346
Show file tree
Hide file tree
Showing 4 changed files with 94 additions and 160 deletions.
2 changes: 1 addition & 1 deletion modules/elasticity/ElementMatrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@
*
* Theory:
*
* a(𝐮,𝐯) = ∫∫ [σ(𝐮):ε(𝐯)dΩ with 𝐮 = (𝑢𝑥,𝑢𝑦) and 𝐯 = (𝑣𝑥,𝑣𝑦)
* a(𝐮,𝐯) = ∫∫ σ(𝐮):ε(𝐯)dΩ with 𝐮 = (𝑢𝑥,𝑢𝑦) and 𝐯 = (𝑣𝑥,𝑣𝑦)
* σ(𝐮) is stress tensor with σᵢⱼ = λδᵢⱼεₖₖ + 2μεᵢⱼ
* ε(𝐯) is strain tensor with εᵢⱼ = 0.5 (∂𝑣ᵢ/∂xⱼ + ∂𝑣ⱼ/∂xᵢ)
*
Expand Down
31 changes: 27 additions & 4 deletions modules/elastodynamics/ElementMatrix.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,29 @@
/*---------------------------------------------------------------------------*/

/*---------------------------------------------------------------------------*/
/**
* @brief Computes the element matrix for a triangular element (ℙ1 FE).
*
* Theory:
*
* a(𝐮,𝐯) = ∫∫ [(∂²𝐮/∂𝑡²).(𝐯)]dΩ + ∫∫ [σ(𝐮):ε(𝐯)]dΩ
*
* with trial func 𝐮 = (𝑢𝑥,𝑢𝑦) and test func 𝐯 = (𝑣𝑥,𝑣𝑦),
* σ(𝐮) is stress tensor with σᵢⱼ = λδᵢⱼεₖₖ + 2μεᵢⱼ
* ε(𝐯) is strain tensor with εᵢⱼ = 0.5 (∂𝑣ᵢ/∂xⱼ + ∂𝑣ⱼ/∂xᵢ)
*
* the bilinear integral after Newmark-Beta and damping terms expands to:
*
* a(𝐮,𝐯) = ∫∫ (c₀)(𝐮.𝐯)
* + ∫∫ (c₁)(∂𝑢𝑦/∂𝑦 ∂𝑣𝑥/∂𝑥 + ∂𝑢𝑥/∂𝑥 ∂𝑣𝑦/∂𝑦)
* + ∫∫ (c₁+2c₂)(∂𝑢𝑥/∂𝑥 ∂𝑣𝑥/∂𝑥 + ∂𝑢𝑦/∂𝑦 ∂𝑣𝑦/∂𝑦)
* + ∫∫ (c₂)(∂𝑢𝑦/∂𝑥 + ∂𝑢𝑥/∂𝑦)(∂𝑣𝑥/∂𝑦 + ∂𝑣𝑦/∂𝑥)
*
* with c₀ = (ρ)/(β δ𝑡²) + (ηₘ ρ γ)/(β δ𝑡)
* c₁ = λ + (λ ηₖ γ)/(β δ𝑡)
* c₂ = 2μ + (2μ ηₖ γ)/(β δ𝑡)
*
*/
/*---------------------------------------------------------------------------*/

FixedMatrix<6, 6> FemModule::
Expand All @@ -29,10 +52,10 @@ _computeElementMatrixTRIA3(Cell cell)
FixedMatrix<1, 6> dyUy = { 0., dyu[0], 0., dyu[1], 0., dyu[2] };
IdentityMatrix<6> I6;

FixedMatrix<6, 6> int_Omega_i = c0 / 12. * ((Uy ^ Uy) + (Ux ^ Ux) + I6) * area +
c1 * ((dyUy ^ dxUx) + (dxUx ^ dyUy)) * area +
(c2 + c1) * ((dxUx ^ dxUx) + (dyUy ^ dyUy)) * area +
(c2 / 2) * ((dxUy + dyUx) ^ (dyUx + dxUy)) * area;
FixedMatrix<6, 6> int_Omega_i = (c0 / 12.) * ((Uy ^ Uy) + (Ux ^ Ux) + I6) * area +
(c1) * ((dyUy ^ dxUx) + (dxUx ^ dyUy)) * area +
(2*c2 + c1) * ((dxUx ^ dxUx) + (dyUy ^ dyUy)) * area +
(c2) * ((dxUy + dyUx) ^ (dyUx + dxUy)) * area;

return int_Omega_i;
}
149 changes: 35 additions & 114 deletions modules/elastodynamics/FemModule.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -121,8 +121,6 @@ _getParameters()
if( options()->lambda.isPresent())
lambda = options()->lambda;

mu2 = mu*2; // lame parameter mu * 2

//----- time discretization Newmark-Beta or Generalized-alpha -----//
if (options()->timeDiscretization == "Newmark-beta") {

Expand All @@ -133,7 +131,7 @@ _getParameters()

c0 = rho/(beta*dt*dt) + etam*rho*gamma/beta/dt ;
c1 = lambda + lambda*etak*gamma/beta/dt ;
c2 = 2.*mu + 2.*mu*etak*gamma/beta/dt ;
c2 = mu + mu*etak*gamma/beta/dt ;
c3 = rho/beta/dt - etam*rho*(1-gamma/beta) ;
c4 = rho*( (1.-2.*beta)/2./beta - etam*dt*(1.-gamma/2/beta)) ;
c5 = -lambda*etak*gamma/beta/dt ;
Expand All @@ -154,7 +152,7 @@ _getParameters()

c0 = rho*(1.-alpm)/(beta*dt*dt) + etam*rho*gamma*(1-alpf)/beta/dt ;
c1 = lambda*(1.-alpf) + lambda*etak*gamma*(1.-alpf)/beta/dt ;
c2 = 2.*mu*(1.-alpf) + 2.*mu*etak*gamma*(1.-alpf)/beta/dt ;
c2 = mu*(1.-alpf) + mu*etak*gamma*(1.-alpf)/beta/dt ;
c3 = rho*(1.-alpm)/beta/dt - etam*rho*(1-gamma*(1-alpf)/beta) ;
c4 = rho*( (1.-alpm)*(1.-2.*beta)/2./beta - alpm - etam*dt*(1.-alpf)*(1.-gamma/2/beta)) ;
c5 = lambda*alpf - lambda*etak*gamma*(1.-alpf)/beta/dt ;
Expand Down Expand Up @@ -276,65 +274,8 @@ _assembleLinearOperator()
Real3 m1 = m_node_coord[cell.nodeId(1)];
Real3 m2 = m_node_coord[cell.nodeId(2)];

/*
Real f0 = m_U[cell.nodeId(0)].x;
Real f1 = m_U[cell.nodeId(1)].x;
Real f2 = m_U[cell.nodeId(2)].x;
Real detA = ( m0.x*(m1.y - m2.y) - m0.y*(m1.x - m2.x) + (m1.x*m2.y - m2.x*m1.y) );
Real2 DXU1;
DXU1.x = ( m0.x*(f1 - f2) - f0*(m1.x - m2.x) + (f2*m1.x - f1*m2.x) ) / (2.*area);/// detA;
DXU1.y = ( f0*(m1.y - m2.y) - m0.y*(f1 - f2) + (f1*m2.y - f2*m1.y) ) / (2.*area);/// detA;
f0 = m_U[cell.nodeId(0)].y;
f1 = m_U[cell.nodeId(1)].y;
f2 = m_U[cell.nodeId(2)].y;
Real2 DXU2;
DXU2.x = ( m0.x*(f1 - f2) - f0*(m1.x - m2.x) + (f2*m1.x - f1*m2.x) ) / (2.*area);/// detA;
DXU2.y = ( f0*(m1.y - m2.y) - m0.y*(f1 - f2) + (f1*m2.y - f2*m1.y) ) / (2.*area);/// detA;
f0 = m_V[cell.nodeId(0)].x;
f1 = m_V[cell.nodeId(1)].x;
f2 = m_V[cell.nodeId(2)].x;
Real2 DXV1;
DXV1.x = ( m0.x*(f1 - f2) - f0*(m1.x - m2.x) + (f2*m1.x - f1*m2.x) ) / (2.*area);/// detA;
DXV1.y = ( f0*(m1.y - m2.y) - m0.y*(f1 - f2) + (f1*m2.y - f2*m1.y) ) / (2.*area);/// detA;
f0 = m_V[cell.nodeId(0)].y;
f1 = m_V[cell.nodeId(1)].y;
f2 = m_V[cell.nodeId(2)].y;
Real2 DXV2;
DXV2.x = ( m0.x*(f1 - f2) - f0*(m1.x - m2.x) + (f2*m1.x - f1*m2.x) ) / (2.*area);/// detA;
DXV2.y = ( f0*(m1.y - m2.y) - m0.y*(f1 - f2) + (f1*m2.y - f2*m1.y) ) / (2.*area);/// detA;
f0 = m_A[cell.nodeId(0)].x;
f1 = m_A[cell.nodeId(1)].x;
f2 = m_A[cell.nodeId(2)].x;
Real2 DXA1;
DXA1.x = ( m0.x*(f1 - f2) - f0*(m1.x - m2.x) + (f2*m1.x - f1*m2.x) ) / (2.*area);/// detA;
DXA1.y = ( f0*(m1.y - m2.y) - m0.y*(f1 - f2) + (f1*m2.y - f2*m1.y) ) / (2.*area);/// detA;
f0 = m_A[cell.nodeId(0)].y;
f1 = m_A[cell.nodeId(1)].y;
f2 = m_A[cell.nodeId(2)].y;
Real2 DXA2;
DXA2.x = ( m0.x*(f1 - f2) - f0*(m1.x - m2.x) + (f2*m1.x - f1*m2.x) ) / (2.*area);/// detA;
DXA2.y = ( f0*(m1.y - m2.y) - m0.y*(f1 - f2) + (f1*m2.y - f2*m1.y) ) / (2.*area);/// detA;
*/

Real2 DXU1, DXU2, DXV1, DXV2, DXA1, DXA2;

// Real2 Ctriangle;
// Ctriangle.x = (1/3.)* (m0.x + m1.x+m2.x);
// Ctriangle.y = (1/3.)* (m0.y + m1.y+m2.y);

// to construct dx(v) we use d(Phi0)/dx , d(Phi1)/dx , d(Phi1)/dx
// here Phi_i are the basis functions at three nodes i=1:3
Real2 dPhi0(m1.y - m2.y, m2.x - m1.x);
Expand All @@ -343,25 +284,24 @@ _assembleLinearOperator()

FixedMatrix<1, 3> DYV;

DYV(0,0) = dPhi0.y /(2.* area);
DYV(0,1) = dPhi1.y /(2.* area);
DYV(0,2) = dPhi2.y /(2.* area);
DYV(0, 0) = dPhi0.y / (2. * area);
DYV(0, 1) = dPhi1.y / (2. * area);
DYV(0, 2) = dPhi2.y / (2. * area);

FixedMatrix<1, 3> DXV;

DXV(0,0) = dPhi0.x /(2.* area);
DXV(0,1) = dPhi1.x /(2.* area);
DXV(0,2) = dPhi2.x /(2.* area);

DXV(0, 0) = dPhi0.x / (2. * area);
DXV(0, 1) = dPhi1.x / (2. * area);
DXV(0, 2) = dPhi2.x / (2. * area);

// to construct dx(u_n) we use d(Phi0)/dx , d(Phi1)/dx , d(Phi1)/dx
// d(u_n)/dx = \sum_{1=1}^{3} { (u_n)_i* (d(Phi_i)/dx) }
Real f0 = m_U[cell.nodeId(0)].x;
Real f1 = m_U[cell.nodeId(1)].x;
Real f2 = m_U[cell.nodeId(2)].x;

DXU1.x = DXV(0,0) * f0 + DXV(0,1) * f1 + DXV(0,2) * f2;
DXU1.y = DYV(0,0) * f0 + DYV(0,1) * f1 + DYV(0,2) * f2;
DXU1.x = DXV(0, 0) * f0 + DXV(0, 1) * f1 + DXV(0, 2) * f2;
DXU1.y = DYV(0, 0) * f0 + DYV(0, 1) * f1 + DYV(0, 2) * f2;

Real Uold1 = f0 + f1 + f2;

Expand All @@ -371,67 +311,52 @@ _assembleLinearOperator()

Real Uold2 = f0 + f1 + f2;

DXU2.x = DXV(0,0) * f0 + DXV(0,1) * f1 + DXV(0,2) * f2;
DXU2.y = DYV(0,0) * f0 + DYV(0,1) * f1 + DYV(0,2) * f2;
DXU2.x = DXV(0, 0) * f0 + DXV(0, 1) * f1 + DXV(0, 2) * f2;
DXU2.y = DYV(0, 0) * f0 + DYV(0, 1) * f1 + DYV(0, 2) * f2;

f0 = m_V[cell.nodeId(0)].x;
f1 = m_V[cell.nodeId(1)].x;
f2 = m_V[cell.nodeId(2)].x;

Real Vold1 = f0 + f1 + f2;

DXV1.x = DXV(0,0) * f0 + DXV(0,1) * f1 + DXV(0,2) * f2;
DXV1.y = DYV(0,0) * f0 + DYV(0,1) * f1 + DYV(0,2) * f2;
DXV1.x = DXV(0, 0) * f0 + DXV(0, 1) * f1 + DXV(0, 2) * f2;
DXV1.y = DYV(0, 0) * f0 + DYV(0, 1) * f1 + DYV(0, 2) * f2;

f0 = m_V[cell.nodeId(0)].y;
f1 = m_V[cell.nodeId(1)].y;
f2 = m_V[cell.nodeId(2)].y;

Real Vold2 = f0 + f1 + f2;

DXV2.x = DXV(0,0) * f0 + DXV(0,1) * f1 + DXV(0,2) * f2;
DXV2.y = DYV(0,0) * f0 + DYV(0,1) * f1 + DYV(0,2) * f2;
DXV2.x = DXV(0, 0) * f0 + DXV(0, 1) * f1 + DXV(0, 2) * f2;
DXV2.y = DYV(0, 0) * f0 + DYV(0, 1) * f1 + DYV(0, 2) * f2;

f0 = m_A[cell.nodeId(0)].x;
f1 = m_A[cell.nodeId(1)].x;
f2 = m_A[cell.nodeId(2)].x;

Real Aold1 = f0 + f1 + f2;

DXA1.x = DXV(0,0) * f0 + DXV(0,1) * f1 + DXV(0,2) * f2;
DXA1.y = DYV(0,0) * f0 + DYV(0,1) * f1 + DYV(0,2) * f2;
DXA1.x = DXV(0, 0) * f0 + DXV(0, 1) * f1 + DXV(0, 2) * f2;
DXA1.y = DYV(0, 0) * f0 + DYV(0, 1) * f1 + DYV(0, 2) * f2;

f0 = m_A[cell.nodeId(0)].y;
f1 = m_A[cell.nodeId(1)].y;
f2 = m_A[cell.nodeId(2)].y;

Real Aold2 = f0 + f1 + f2;

DXA2.x = DXV(0,0) * f0 + DXV(0,1) * f1 + DXV(0,2) * f2;
DXA2.y = DYV(0,0) * f0 + DYV(0,1) * f1 + DYV(0,2) * f2;


/*
$$
\int_{\Omega}(
(U \cdot v) c_0
+ (V \cdot v) c_3
+ (A \cdot v) c_4
- (\nabla \cdot U \nabla \cdot v) c_5
- (\varepsilon(U) : \varepsilon(v) ) c_6
+ (\nabla \cdot V \nabla \cdot v) c_7
+ (\varepsilon(V) : \varepsilon(v) ) c_9
+ (\nabla \cdot A \nabla \cdot v) c_8
+ (\varepsilon(A) : \varepsilon(v) ) c_{10}
)
$$
*/

// Info: We could also use the following logic
// cell.node(i).isOwn();
// cell.node(i);
// Node node = cell.node(i);
// Then we can loop 1:3
DXA2.x = DXV(0, 0) * f0 + DXV(0, 1) * f1 + DXV(0, 2) * f2;
DXA2.y = DYV(0, 0) * f0 + DYV(0, 1) * f1 + DYV(0, 2) * f2;

// the linear terms expands to:
//
// b(0,𝐯) = ∫∫ (c₀)(𝐮ₙ.𝐯) + ∫∫ (c₃)(ụₙ.𝐯) + ∫∫ (c₄)(ṳₙ.𝐯)
// - ∫∫ (c₅)(∇𝐮ₙ.∇𝐯) - ∫∫ (c₆)(ε(𝐮ₙ):ε(𝐯))
// + ∫∫ (c₇)(∇ụₙ.∇𝐯) + ∫∫ (c₈)(ε(ụₙ):ε(𝐯))
// + ∫∫ (c₉)(∇ṳₙ.∇𝐯) + ∫∫ (c₁₀)(ε(ṳₙ):ε(𝐯))

int i = 0;
for (Node node : cell.nodes()) {
if (node.isOwn()) {
Expand Down Expand Up @@ -499,28 +424,24 @@ _assembleLinearOperator()
Real length = ArcaneFemFunctions::MeshOperation::computeLengthEdge2(face, m_node_coord);
for (Node node : iface->nodes()) {
if (node.isOwn()) {
DoFLocalId dof_id1 = node_dof.dofId(node, 0);
rhs_values[dof_id1] += trac.x * length / 2.;
}
if (node.isOwn()) {
DoFLocalId dof_id2 = node_dof.dofId(node, 1);
rhs_values[dof_id2] += trac.y * length / 2.;
rhs_values[node_dof.dofId(node, 0)] += trac.x * length / 2.;
rhs_values[node_dof.dofId(node, 1)] += trac.y * length / 2.;
}
}
}
continue;
}
else {
const UniqueArray<String> t_string = bs->t();
Real3 t;
Real3 trac;

info() << "[ArcaneFem-Info] Applying Traction " << t_string;
info() << "[ArcaneFem-Info] Traction surface '" << bs->surface().name() << "'";

for (Int32 i = 0; i < t_string.size(); ++i) {
t[i] = 0.0;
trac[i] = 0.0;
if (t_string[i] != "NULL") {
t[i] = std::stod(t_string[i].localstr());
trac[i] = std::stod(t_string[i].localstr());
}
}

Expand All @@ -531,8 +452,8 @@ _assembleLinearOperator()
Real length = ArcaneFemFunctions::MeshOperation::computeLengthEdge2(face, m_node_coord);
for (Node node : iface->nodes()) {
if (node.isOwn()) {
rhs_values[node_dof.dofId(node, 0)] += t[0] * length / 2.;
rhs_values[node_dof.dofId(node, 1)] += t[1] * length / 2.;
rhs_values[node_dof.dofId(node, 0)] += trac[0] * length / 2.;
rhs_values[node_dof.dofId(node, 1)] += trac[1] * length / 2.;
}
}
}
Expand Down
72 changes: 31 additions & 41 deletions modules/elastodynamics/FemModule.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,50 +60,40 @@ class FemModule
delete t.case_table;
}

public:

//! Method called at each iteration
void compute() override;

//! Method called at the beginning of the simulation
void startInit() override;

VersionInfo versionInfo() const override
{
return VersionInfo(1, 0, 0);
}
void startInit() override; //! Method called at the beginning of the simulation
void compute() override; //! Method called at each iteration
VersionInfo versionInfo() const override { return VersionInfo(1, 0, 0); }

private:

Real t; // time variable
Real dt; // time step
Real tmax; // max time

Real etam; // time discretization param etam
Real etak; // time discretization param etak
Real alpm; // time discretization param alpm
Real alpf; // time discretization param alpf
Real beta; // time discretization param beta
Real gamma; // time discretization param gamma

Real E; // Youngs modulus
Real nu; // Poissons ratio
Real rho; // Density
Real mu; // Lame parameter mu
Real mu2; // Lame parameter mu * 2
Real lambda; // Lame parameter lambda

Real c0; // constant
Real c1; // constant
Real c2; // constant
Real c3; // constant
Real c4; // constant
Real c5; // constant
Real c6; // constant
Real c7; // constant
Real c8; // constant
Real c9; // constant
Real c10; // constant
Real t; // time variable 𝑡
Real dt; // time step δ𝑡
Real tmax; // max time 𝑡ₘₐₓ
Real alpm; // time discretization param αᵐ
Real alpf; // time discretization param αᶠ
Real beta; // time discretization param β
Real gamma; // time discretization param γ

Real etam; // damping parameter ηₘ
Real etak; // damping parameter ηₖ

Real E; // Youngs modulus 𝐸
Real nu; // Poissons ratio ν
Real rho; // Density ρ
Real mu; // Lame parameter μ
Real lambda; // Lame parameter λ

Real c0; // constant c₀
Real c1; // constant c₁
Real c2; // constant c₂
Real c3; // constant c₃
Real c4; // constant c₄
Real c5; // constant c₅
Real c6; // constant c₆
Real c7; // constant c₇
Real c8; // constant c₈
Real c9; // constant c₉
Real c10; // constant c₁₀

DoFLinearSystem m_linear_system;
FemDoFsOnNodes m_dofs_on_nodes;
Expand Down

0 comments on commit 51c9346

Please sign in to comment.