3D soild structure eigen problem

Léonhard YU

unread,

Mar 28, 2022, 11:09:28 PM3/28/22

to deal.II User Group

Dear Friends,

I am constructing a code to solve the eigen problem of a solid structure,.

I use the step-36 as a sample, but I don't know how to write a new one for 3D solid structure with tensil modulus E = 2.0e11, poisson ratio = 0.3 and density = 7.85e3.

If you have any other samples or ideas, please help me.

Thanks a lot,

best wishes

Wolfgang Bangerth

unread,

Mar 29, 2022, 10:20:30 AM3/29/22

to dea...@googlegroups.com

On 3/28/22 21:09, Léonhard YU wrote:
> I am constructing a code to solve the eigen problem of a solid structure,.
> I use the step-36 as a sample, but I don't know how to write a new one for 3D
> solid structure with tensil modulus E = 2.0e11, poisson ratio = 0.3 and
> density = 7.85e3.
> If you have any other samples or ideas, please help me.

Leonhard:
what have you tried already? Have you looked at the elasticity tutorial
programs and how they construct (i) finite element spaces and (ii) the
stiffness matrix? If you understand that, you will also understand how to
solve the eigenproblem.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

Léonhard YU

unread,

Mar 29, 2022, 9:26:45 PM3/29/22

to deal.II User Group

Thanks for your reply, professor.

I have constructed the stiffness matrix as step-18 does, but the difficulty is the mass matrix. I use these to construct mass matrix:

1) cell_mass_matrix(i, j) += rho * fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point) * fe_values.JxW(q_point);

2) cell_mass_matrix(i, j) +=
(mass_phi_i * mass_phi_j) *
fe_values.JxW(q_point) * rho;

where mass_phi_i:

const Tensor<2, dim>
mass_phi_i =get_mass_matrix(fe_values, i, q_point);

(using this function)

template <int dim>

inline Tensor<2, dim> get_mass_matrix(const FEValues<dim> &fe_values,
const unsigned int shape_func_i,
const unsigned int shape_func_j,
const unsigned int q_point)
{
Tensor<2, dim> tmp;
for (unsigned int t = 0; t < dim; ++t)
for (unsigned int k = 0; k < dim; ++k)
{
tmp[t][k] = fe_values.shape_value_component(shape_func_i, q_point, t)*
fe_values.shape_value_component(shape_func_j, q_point, k);
}
return tmp;
}

2) cell_mass_matrix(i, j) +=
(mass_phi_i * mass_phi_j) *
fe_values.JxW(q_point) * rho;

where mass_phi_i:

const SymmetricTensor <2, dim>
mass_phi_i =get_mass_matrix(fe_values, i, q_point);

(using this function)

template <int dim>
inline SymmetricTensor<2, dim>
get_mass_matrix (const FEValues<dim> &fe_values,
const unsigned int shape_func,
const unsigned int q_point)
{
SymmetricTensor<2, dim> tmp;
for (unsigned int i = 0; i < dim; ++i)
tmp[i][i] = fe_values.shape_value_component(shape_func, q_point, i) *
fe_values.shape_value_component(shape_func, q_point, i);

for (unsigned int i = 0; i < dim; ++i)
for (unsigned int j = i + 1; j < dim; ++j)
tmp[i][j] =
fe_values.shape_value_component(shape_func, q_point, i) *
fe_values.shape_value_component(shape_func, q_point, j);
return tmp;
}

Either of them is not right.

Wolfgang Bangerth

unread,

Mar 30, 2022, 4:10:24 PM3/30/22

to dea...@googlegroups.com

On 3/29/22 19:26, Léonhard YU wrote:
> 1) cell_mass_matrix(i, j) += rho * fe_values.shape_value(i, q_point) *
> fe_values.shape_value(j, q_point) * fe_values.JxW(q_point);

Yes, this is not correct unless you prefix it with
if (fe.system_to_component_index(i).first ==
fe.system_to_component_index(j).first)

> 2) cell_mass_matrix(i, j) +=
> (mass_phi_i * mass_phi_j) *
> fe_values.JxW(q_point) * rho;
>
> where mass_phi_i:
> const Tensor<2, dim>
> mass_phi_i =get_mass_matrix(fe_values, i, q_point);
>
> (using this function)
> template <int dim>
> inline Tensor<2, dim> get_mass_matrix(const FEValues<dim> &fe_values,
> const unsigned int shape_func_i,
> const unsigned int shape_func_j,
> const unsigned int q_point)
> {
> Tensor<2, dim> tmp;
> for (unsigned int t = 0; t < dim; ++t)
> for (unsigned int k = 0; k < dim; ++k)
> {
> tmp[t][k] = fe_values.shape_value_component(shape_func_i,
> q_point, t)*
> fe_values.shape_value_component(shape_func_j,
> q_point, k);
> }
> return tmp;
> }

I don't understand this code. You say

const Tensor<2, dim> mass_phi_i
= get_mass_matrix(fe_values, i, q_point);

but the definition of get_mass_matrix() takes four arguments.

> 2) cell_mass_matrix(i, j) +=
> (mass_phi_i * mass_phi_j) *
> fe_values.JxW(q_point) * rho;
>
> where mass_phi_i:
> const SymmetricTensor <2, dim>
> mass_phi_i =get_mass_matrix(fe_values, i, q_point);
>
> (using this function)
> template <int dim>
> inline SymmetricTensor<2, dim>
> get_mass_matrix (const FEValues<dim> &fe_values,
> const unsigned int shape_func,
> const unsigned int q_point)
> {
> SymmetricTensor<2, dim> tmp;
> for (unsigned int i = 0; i < dim; ++i)
> tmp[i][i] = fe_values.shape_value_component(shape_func, q_point, i) *
> fe_values.shape_value_component(shape_func,
> q_point, i);
>
> for (unsigned int i = 0; i < dim; ++i)
> for (unsigned int j = i + 1; j < dim; ++j)
> tmp[i][j] =
> fe_values.shape_value_component(shape_func, q_point, i) *
> fe_values.shape_value_component(shape_func, q_point, j);
> return tmp;
> }

This suffers from the same problem.

I think you want to read through the module on vector-valued problems to
understand the general philosophy about dealing with vector problems.

Léonhard YU

unread,

Apr 6, 2022, 7:46:29 AM4/6/22

to deal.II User Group

Dear Prof. Bangerth,

Thanks for your great help.

I made a mistake on copying the definition of function at 2), but it is not right either.

I rewrote the codes.

The stiffness matrix is constructed according to (Latex):

\left[\mathbf{k}_{m}\right]=\iiint[\mathbf{B}]^{T}[\mathbf{D}][\mathbf{B}] d x d y d z

where:

[\mathbf{D}]=\frac{E(1-v)}{(1+v)(1-2 v)}\left[\begin{array}{ccccccc}
1 & \frac{v}{1-v} & \frac{v}{1-v} & 0 & 0 & 0 \\
\frac{v}{1-v} & 1 & \frac{v}{1-v} & 0 & 0 & 0 \\
\frac{v}{1-v} & \frac{v}{1-v} & 1 & 0 & 0 & 0 \\
0 & 0 & 0 & \frac{1-2 v}{2(1-v)} & 0 & 0 \\
0 & 0 & 0 & 0 & \frac{1-2 v}{2(1-v)} & 0 \\
0 & 0 & 0 & 0 & 0 & \frac{1-2 v}{2(1-v)}
\end{array}\right]

[\mathbf{B}]=[\mathbf{A}][\mathbf{S}]

[\mathbf{A}]=\left[\begin{array}{ccc}
\frac{\partial}{\partial x} & 0 & 0 \\
0 & \frac{\partial}{\partial y} & 0 \\
0 & 0 & \frac{\partial}{\partial z} \\
\frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \\
0 & \frac{\partial}{\partial z} & \frac{\partial}{\partial y} \\
\frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial x}
\end{array}\right]

[\mathbf{S}]=\left[\begin{array}{cccccccccc}
N_{1} & 0 & 0 & N_{2} & 0 & 0 & N_{3} & 0 & 0 & \cdots \\
0 & N_{1} & 0 & 0 & N_{2} & 0 & 0 & N_{3} & 0 & \cdots \\
0 & 0 & N_{1} & 0 & 0 & N_{2} & 0 & 0 & N_{3} & \cdots
\end{array}\right]

and mass matrix:

\left[\mathbf{m}_{m}\right] = \rho \iiint[\mathbf{N}]^{T}[\mathbf{N}] d x d y d z

The number of degree of freedom of 8-node cell is 8 according to the reference book, but in deal.II it is 8*3.

I do not understand that. I guess that each column in matrix S is for one DoF, but I am not sure.

The cell stiffness and mass matrix are:

for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
{
Tensor<1,6> eps_u_i;
Tensor<1,6> eps_u_j;
switch (fe.system_to_component_index(i).first){
case 0:
{
eps_u_i[0] = fe_values.shape_grad_component(i, q_point, 0)[0];
eps_u_i[3] = fe_values.shape_grad_component(i, q_point, 0)[1];
eps_u_i[5] = fe_values.shape_grad_component(i, q_point, 0)[2];
break;
}
case 1:
{
eps_u_i[1] = fe_values.shape_grad_component(i, q_point, 1)[1];
eps_u_i[3] = fe_values.shape_grad_component(i, q_point, 1)[0];
eps_u_i[4] = fe_values.shape_grad_component(i, q_point, 1)[2];
break;
}
case 2:
{
eps_u_i[2] = fe_values.shape_grad_component(i, q_point, 2)[2];
eps_u_i[4] = fe_values.shape_grad_component(i, q_point, 2)[1];
eps_u_i[5] = fe_values.shape_grad_component(i, q_point, 2)[0];
break;
}
default:
{
std::cout << "wrong here" << std::endl;
break;
}
}

switch (fe.system_to_component_index(j).first){
case 0:
{
eps_u_j[0] = fe_values.shape_grad_component(j, q_point, 0)[0];
eps_u_j[3] = fe_values.shape_grad_component(j, q_point, 0)[1];
eps_u_j[5] = fe_values.shape_grad_component(j, q_point, 0)[2];
break;
}
case 1:
{
eps_u_j[1] = fe_values.shape_grad_component(j, q_point, 1)[1];
eps_u_j[3] = fe_values.shape_grad_component(j, q_point, 1)[0];
eps_u_j[4] = fe_values.shape_grad_component(j, q_point, 1)[2];
break;
}
case 2:
{
eps_u_j[2] = fe_values.shape_grad_component(j, q_point, 2)[2];
eps_u_j[4] = fe_values.shape_grad_component(j, q_point, 2)[1];
eps_u_j[5] = fe_values.shape_grad_component(j, q_point, 2)[0];
break;
}
default:
{
std::cout << "wrong here" << std::endl;
break;
}
}
cell_stiffness_matrix(i, j) *= eps_u_i * stress_strain_D * eps_u_j * fe_values.JxW(q_point);

if (fe.system_to_component_index(i).first ==
fe.system_to_component_index(j).first)
{

cell_mass_matrix(i, j) +=
rho *
fe_values.shape_value(i, q_point) *
fe_values.shape_value(j, q_point) *
fe_values.JxW(q_point);

}
}

cell->get_dof_indices(local_dof_indices);
constraints.distribute_local_to_global(
cell_stiffness_matrix, local_dof_indices, stiffness_matrix);
constraints.distribute_local_to_global(cell_mass_matrix,
local_dof_indices,
mass_matrix);
}

Now it reports convergence only after 1 step, the eigenvalues are all drop in the interval of spurious eigenvalues.

Number of active cells: 4875
Number of degrees of freedom: 19008
Spurious eigenvalues are all in the interval [1.11783e-05,1.11783e-05]
Solver converged in 1 iterations.

Eigenvalue 0 : 1.11783e-05
Eigenvalue 1 : 1.11783e-05
Eigenvalue 2 : 1.11783e-05
Eigenvalue 3 : 1.11783e-05
Eigenvalue 4 : 1.11783e-05
Eigenvalue 5 : 1.11783e-05
Eigenvalue 6 : 1.11783e-05
Eigenvalue 7 : 1.11783e-05
Eigenvalue 8 : 1.11783e-05
Eigenvalue 9 : 1.11783e-05
Eigenvalue 10 : 1.11783e-05
Eigenvalue 11 : 1.11783e-05
Eigenvalue 12 : 1.11783e-05
Eigenvalue 13 : 1.11783e-05
Eigenvalue 14 : 1.11783e-05
Eigenvalue 15 : 1.11783e-05
Eigenvalue 16 : 1.11783e-05
Eigenvalue 17 : 1.11783e-05
Eigenvalue 18 : 1.11783e-05
Eigenvalue 19 : 1.11783e-05
[100%] Built target run

I use the codes

for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
if (constraints.is_constrained(i))
{
stiffness_matrix.set(i, i, 1.11783e-05);
mass_matrix.set(i, i, 1);
}

as step-36 to deal with it, but errors are:

[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Not for unassembled matrix
[0]PETSC ERROR: See https://petsc.org/release/faq/ for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.16.3, Jan 05, 2022
[0]PETSC ERROR: ./Test4 on a arch-linux-c-debug named eineschraube-virtual-machine by eineschraube Wed Apr 6 17:07:57 2022
[0]PETSC ERROR: Configure options --with-mpidir=/usr/lib --download-fblaslapack
[0]PETSC ERROR: #1 MatGetOrdering() at /home/eineschraube/Code/petsc-3.16.3/src/mat/order/sorder.c:177
[0]PETSC ERROR: #2 PCSetUp_LU() at /home/eineschraube/Code/petsc-3.16.3/src/ksp/pc/impls/factor/lu/lu.c:88
[0]PETSC ERROR: #3 PCSetUp() at /home/eineschraube/Code/petsc-3.16.3/src/ksp/pc/interface/precon.c:1016
[0]PETSC ERROR: #4 KSPSetUp() at /home/eineschraube/Code/petsc-3.16.3/src/ksp/ksp/interface/itfunc.c:408
[0]PETSC ERROR: #5 STSetUp_Shift() at /home/eineschraube/Code/slepc-3.16.2/src/sys/classes/st/impls/shift/shift.c:107
[0]PETSC ERROR: #6 STSetUp() at /home/eineschraube/Code/slepc-3.16.2/src/sys/classes/st/interface/stsolve.c:582
[0]PETSC ERROR: #7 EPSSetUp() at /home/eineschraube/Code/slepc-3.16.2/src/eps/interface/epssetup.c:350
[0]PETSC ERROR: #8 EPSSolve() at /home/eineschraube/Code/slepc-3.16.2/src/eps/interface/epssolve.c:136
terminate called after throwing an instance of 'dealii::SLEPcWrappers::SolverBase::ExcSLEPcError'
what():
--------------------------------------------------------
An error occurred in line <185> of file </home/eineschraube/Code/dealii/dealii-9.4.0pre/source/lac/slepc_solver.cc> in function
void dealii::SLEPcWrappers::SolverBase::solve(unsigned int, unsigned int*)
The violated condition was:
ierr == 0
Additional information:
An error with error number 73 occurred while calling a SLEPc function

Stacktrace:
-----------
#0 /home/eineschraube/Code/dealii/dealii-9.4.0pre/build/lib/libdeal_II.g.so.9.4.0-pre: dealii::SLEPcWrappers::SolverBase::solve(unsigned int, unsigned int*)
#1 ./Test4: void dealii::SLEPcWrappers::SolverBase::solve<dealii::PETScWrappers::MPI::Vector>(dealii::PETScWrappers::MatrixBase const&, dealii::PETScWrappers::MatrixBase const&, std::vector<double, std::allocator<double> >&, std::vector<dealii::PETScWrappers::MPI::Vector, std::allocator<dealii::PETScWrappers::MPI::Vector> >&, unsigned int)
#2 ./Test4: Test4::BeamEigenProblem<3>::solve()
#3 ./Test4: Test4::BeamEigenProblem<3>::run()
#4 ./Test4: main

I do know how to handle this question.

Sincere thank for your reply and help.

Best wishes,

YU

Wolfgang Bangerth

unread,

Apr 6, 2022, 1:38:12 PM4/6/22

to dea...@googlegroups.com

On 4/6/22 05:46, Léonhard YU wrote:
>
> The number of degree of freedom of 8-node cell is 8 according to the
> reference book, but in deal.II it is 8*3.
> I do not understand that. I guess that each column in matrix S is for
> one DoF, but I am not sure.

I will point you again at the vector-valued documentation module that
has many examples. In your case, you have x-, y-, and z-displacements at
each vertex of the cell, so you have 8*3 degrees of freedom on each cell.

> The cell stiffness and mass matrix are:
>
> for (unsigned int i = 0; i < dofs_per_cell; ++i)
> for (unsigned int j = 0; j < dofs_per_cell; ++j)
> for (unsigned int q_point = 0; q_point < n_q_points;
> ++q_point)
> {
> Tensor<1,6> eps_u_i;

This suggests that you are using Voigt notation. This is common in the
engineering literature, but not typically the perspective we take in
(the more mathematically oriented) deal.II. Rather, we tend to use
notation that emphasizes that eps_u_i is a rank-2 symmetric tensor, and
so use SymmetricTensor<2,3> instead.

It is not inherently wrong to do as you do but you will find more
examples in the deal.II tutorials if you use tensor notation.

> cell_stiffness_matrix(i, j) *= eps_u_i * stress_strain_D
> * eps_u_j * fe_values.JxW(q_point);

I am certain that you want to use

> cell_stiffness_matrix(i, j) += eps_u_i * stress_strain_D
> * eps_u_j * fe_values.JxW(q_point);

instead (note the += instead of *=).

> if (fe.system_to_component_index(i).first ==
> fe.system_to_component_index(j).first)
> {
> cell_mass_matrix(i, j) +=
> rho *
> fe_values.shape_value(i, q_point) *
> fe_values.shape_value(j, q_point) *
> fe_values.JxW(q_point);
> }
> }

This looks correct.

> Now it reports convergence only after 1 step, the eigenvalues are all
> drop in the interval of spurious eigenvalues.

Yes. That's because the mistake above means that your stiffness matrix
entries are probably all zero, with the exception of those you set to
something else by hand.

Reply all

Reply to author

Forward