Hi, Tobi:
Thank you so much for your suggestions. I have followed your direction and
checked my code again and again for all day long, but still
could not resolve this problem. If you (or someone else) wouldn't mind,
could you please take a look at my code attached below?
It is a really short code, of which most parts are exactly the same with
step-36, except some details I changed for learning purposes (I changed
this problem to a simple dynamical elasticity problem on a rectangular
domain). I know this might cause some unconveniences, but I really
need somesone's help. Since nobody around me uses dealii, this is the only
place I can turn to.
Thanks in advance for any help!
Best regards,
David.
On Monday, May 16, 2016 at 1:03:30 AM UTC+8, Tobi Young wrote:
>
> > Could anyone please take a look, it says:
> >
> > [0]PETSC ERROR: Object is in wrong state
> > [0]PETSC ERROR: Matrix is missing diagonal entry 8
> >
>
> >The above tells us that the matrix is singular (missing a diagonal
> entry).
>
>
>I can recommend you check first, if your matrix should be singular -
> >ie. if you are assembling him correctly; and second, check that the
> >SLEPc solver you are using is capable of solving a singular matrix
> >(not all are). You can get some information from the solver pages of
> >the SLEPc manual.
> >Hope that helps a little.
>
--
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* ---------------------------------------------------------------------
*
* Copyright (C) 2000 - 2015 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* ---------------------------------------------------------------------
*/
#include <deal.II/base/logstream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/function_parser.h>
#include <deal.II/lac/petsc_sparse_matrix.h>
#include <deal.II/lac/petsc_parallel_vector.h>
#include <deal.II/lac/slepc_solver.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <fstream>
#include <iostream>
namespace eigenproblem
{
using namespace dealii;
//Define material constants.
double Es = 22407.96E6; //Young's modulus
double rho = 2483.; //Mass density
double nu = 0.2; //Poisson's ratio
double xi = 0.05; //Damping ratio
template <int dim, int fe_degree>
class ModesComputation
{
public:
ModesComputation (const FiniteElement<dim> &finite_element);
void run ();
private:
void make_triangulation ();
void setup_system ();
void assemble_system ();
void solve ();
void output_results () const;
Triangulation<dim> triangulation;
FESystem<dim> fe;
DoFHandler<dim> dof_handler;
ConstraintMatrix constraints;
PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
std::vector<PETScWrappers::MPI::Vector> eigenfunctions;
std::vector<double> eigenvalues;
unsigned int eigen_numbers;
double Elasticity_Tensor(unsigned int i,unsigned int j,
unsigned int k,unsigned int l);
};
template <int dim, int fe_degree>
ModesComputation<dim, fe_degree>::ModesComputation(const FiniteElement<dim> &finite_element)
:
fe (finite_element, dim),
dof_handler (triangulation),
eigen_numbers(8)
{}
template <int dim, int fe_degree>
double ModesComputation<dim, fe_degree>::Elasticity_Tensor(unsigned int i,unsigned int j,
unsigned int k,unsigned int l)
{
double lambda=(Es*nu)/((1.+nu)*(1.-2.*nu)),
mu=Es/(2.*(1.+nu));
return lambda*(i==j)*(k==l) + mu*((i==k)*(j==l) + (i==l)*(j==k));
}
template <int dim, int fe_degree>
void ModesComputation<dim, fe_degree>::make_triangulation()
{
Point<2> point1 =Point<2> (0., 0.);
Point<2> point2 =Point<2> (2., 20.);
GridGenerator::hyper_rectangle (triangulation,
point1,
point2,
true);
triangulation.refine_global(4);
std::string filename = "rectangle.vtk";
std::ofstream out (filename.c_str());
GridOut grid_out;
grid_out.write_vtk (triangulation, out);
std::cout << " written to " << filename
<< std::endl
<< std::endl;
}
template <int dim, int fe_degree>
void ModesComputation<dim, fe_degree>::setup_system()
{
dof_handler.distribute_dofs(fe);
std::cout << " Number of degrees of freedom: "
<< dof_handler.n_dofs ()
<< std::endl;
DoFTools::make_zero_boundary_constraints (dof_handler, 2, constraints);
constraints.close();
stiffness_matrix.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
dof_handler.max_couplings_between_dofs());
mass_matrix.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
dof_handler.max_couplings_between_dofs());
IndexSet eigenfunction_index_set = dof_handler.locally_owned_dofs ();
for (unsigned int i=0; i<eigenfunctions.size (); ++i)
eigenfunctions[i].reinit (eigenfunction_index_set, MPI_COMM_WORLD);
eigenfunctions.resize (eigen_numbers);
eigenvalues.resize (eigenfunctions.size ());
}
template <int dim, int fe_degree>
void ModesComputation<dim, fe_degree>::assemble_system()
{
QGaussLobatto<dim> quadrature_formula(fe_degree + 1);
FEValues<dim> fe_values (fe,
quadrature_formula,
update_values |
update_gradients |
update_JxW_values);
stiffness_matrix=0; mass_matrix=0;
const unsigned int dofs_per_elem = fe.dofs_per_cell;
const unsigned int nodes_per_elem = GeometryInfo<dim>::vertices_per_cell;
const unsigned int num_quad_pts = quadrature_formula.size();
FullMatrix<double> Klocal (dofs_per_elem, dofs_per_elem),
Mlocal (dofs_per_elem, dofs_per_elem);
std::vector<unsigned int> local_dof_indices (dofs_per_elem);
typename DoFHandler<dim>::active_cell_iterator elem = dof_handler.begin_active(), endc = dof_handler.end();
for (;elem!=endc; ++elem)
{
fe_values.reinit(elem);
elem->get_dof_indices (local_dof_indices);
Klocal = 0.;
Mlocal = 0.;
for (unsigned int q=0; q<num_quad_pts; ++q){
for (unsigned int A=0; A<nodes_per_elem; A++) {
for(unsigned int i=0; i<dim; i++){
for (unsigned int B=0; B<nodes_per_elem; B++) {
for(unsigned int k=0; k<dim; k++){
for (unsigned int j = 0; j<dim; j++){
for (unsigned int l = 0; l<dim; l++){
Klocal[dim*A+i][dim*B+k] += fe_values.shape_grad(dim*A+i, q)[j] *
Elasticity_Tensor(i, j, k, l) *
fe_values.shape_grad(dim*B+k, q)[l] *
fe_values.JxW(q);
}
}
}
}
}
}
}
for (unsigned int q=0; q<num_quad_pts; ++q){
for (unsigned int A=0; A<nodes_per_elem; A++) {
for(unsigned int i=0; i<dim; i++){
for (unsigned int B=0; B<nodes_per_elem; B++) {
for(unsigned int k=0; k<dim; k++){
Mlocal[dim*A+i][dim*B+k] += rho * fe_values.shape_value(dim*A+i, q) *
fe_values.shape_value(dim*B+k, q) *
fe_values.JxW(q) * (i == k);
}
}
}
}
}
constraints
.distribute_local_to_global (Klocal,
local_dof_indices,
stiffness_matrix);
constraints
.distribute_local_to_global (Mlocal,
local_dof_indices,
mass_matrix);
}
stiffness_matrix.compress (VectorOperation::add);
mass_matrix.compress (VectorOperation::add);
double min_spurious_eigenvalue = std::numeric_limits<double>::max(),
max_spurious_eigenvalue = -std::numeric_limits<double>::max();
for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
if (constraints.is_constrained(i))
{
const double ev = stiffness_matrix(i,i)/mass_matrix(i,i);
min_spurious_eigenvalue = std::min (min_spurious_eigenvalue, ev);
max_spurious_eigenvalue = std::max (max_spurious_eigenvalue, ev);
}
std::cout << " Spurious eigenvalues are all in the interval "
<< "[" << min_spurious_eigenvalue << "," << max_spurious_eigenvalue << "]"
<< std::endl;
}
template <int dim, int fe_degree>
void ModesComputation<dim, fe_degree>::solve()
{
SolverControl solver_control (dof_handler.n_dofs(), 1e-9);
SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control);
eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL);
eigensolver.set_problem_type (EPS_GHEP);
eigensolver.solve (stiffness_matrix, mass_matrix,
eigenvalues, eigenfunctions,
eigenfunctions.size());
for (unsigned int i=0; i<eigenfunctions.size(); ++i)
eigenfunctions[i] /= eigenfunctions[i].linfty_norm ();
std::cout << solver_control.last_step () << std::endl;
}
template <int dim, int fe_degree>
void ModesComputation<dim, fe_degree>::output_results() const
{
DataOut<dim> data_out;
data_out.attach_dof_handler (dof_handler);
for (unsigned int i=0; i<eigenfunctions.size(); ++i)
data_out.add_data_vector (eigenfunctions[i],
std::string("eigenfunction_") +
Utilities::int_to_string(i));
data_out.build_patches ();
std::ofstream output ("eigenvectors.vtk");
data_out.write_vtk (output);
}
template <int dim, int fe_degree>
void ModesComputation<dim, fe_degree>::run()
{
make_triangulation ();
setup_system ();
assemble_system ();
solve ();
output_results ();
std::cout << std::endl;
for (unsigned int i=0; i<eigenvalues.size(); ++i)
std::cout << " Eigenvalue " << i
<< " : " << eigenvalues[i]
<< std::endl;
}
}
int main (int argc, char **argv)
{
try
{
using namespace eigenproblem;
using namespace dealii;
dealii::deallog.depth_console (0);
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
AssertThrow(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)==1,
ExcMessage("This program can only be run in serial, use ./eigenproblem"));
const unsigned int dim = 2;
const unsigned int fe_degree = 2;
FE_Q<dim> fe_lobatto(QGaussLobatto<1> (fe_degree+1));
ModesComputation<dim, fe_degree> modes_computation(fe_lobatto);
modes_computation.run ();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
##
# CMake script for the step-8 tutorial program:
##
# Set the name of the project and target:
SET(TARGET "eigenproblem")
# Declare all source files the target consists of. Here, this is only
# the one step-X.cc file, but as you expand your project you may wish
# to add other source files as well. If your project becomes much larger,
# you may want to either replace the following statement by something like
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
# or switch altogether to the large project CMakeLists.txt file discussed
# in the "CMake in user projects" page accessible from the "User info"
# page of the documentation.
SET(TARGET_SRC
${TARGET}.cc
)
# Usually, you will not need to modify anything beyond this point...
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8)
FIND_PACKAGE(deal.II 8.3 QUIET
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
)
IF(NOT ${deal.II_FOUND})
MESSAGE(FATAL_ERROR "\n"
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to
cmake\n"
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
)
ENDIF()
DEAL_II_INITIALIZE_CACHED_VARIABLES()
PROJECT(${TARGET})
DEAL_II_INVOKE_AUTOPILOT()
