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
Le jeudi 31 mars 2022 à 04:10:24 UTC+8, Wolfgang Bangerth a écrit :
> 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.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: [email protected]
> www: http://www.math.colostate.edu/~bangerth/
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/6064a79c-016a-4d9f-868e-4b3530ed1a3en%40googlegroups.com.
