Can anyone help me to fix this.Sinced there is no continuous hermite
interpolation element.I used direct formule for stiffness and mass matix.
Here i am attaching the code.I used step-36.
--
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].
For more options, visit https://groups.google.com/d/optout.
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.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_q.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/lac/vector.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/utilities.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/function_parser.h>
//#include <deal.II/base/function.h>
//#include <deal.II/base/logstream.h>
//#include <deal.II/base/tensor_function.h>
//#include <deal.II/numerics/matrix_tools.h>
//#include <deal.II/fe/mapping.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/lac/full_matrix.h>
#include <deal.II/base/index_set.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/petsc_parallel_sparse_matrix.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 <iostream>
#include <fstream>
#include <cmath>
#include <vector>
using namespace dealii;
template<int dim>
class Eigen
{
public:
Eigen();
~Eigen();
void generate_mesh(unsigned int numberOfElements);
void deg_of_freedom();
void assembly();
unsigned int solve ();
void output_results () const;
//void define_boundary_conds();
void run();
//Class objects
Triangulation<dim> triangulation; //mesh
FESystem<dim> fe; //FE element
DoFHandler<dim> dof_handler; //Connectivity matrices
double L;
//std::map<unsigned int,double> boundary_values; //Map of dirichlet boundary conditions
//std::vector<double> nodeLocation; //Vector of the x-coordinate of nodes by global dof number
PETScWrappers::SparseMatrix stiffness_matrix, mass_matrix;
std::vector<PETScWrappers::MPI::Vector> eigenfunctions;
std::vector<double> eigenvalues;
ConstraintMatrix constraints;
};
template<int dim>
Eigen<dim>::Eigen()
:
dof_handler(triangulation),/* fe(degree of the polynomial of the element)*/
fe(FE_Q <dim>(1),2) /*each node 2 components,as it is 1d so 1+1*/
{}
template <int dim>
Eigen<dim>::~Eigen(){
dof_handler.clear();
}
template <int dim>
void Eigen<dim>::generate_mesh(unsigned int numberOfElements)
{
//Define the limits of your domain
L =3 ; //EDIT
double x_min = 0.;
double x_max = L;
Point<dim,double> min(x_min),
max(x_max);
std::vector<unsigned int> meshDimensions (dim,numberOfElements);
GridGenerator::subdivided_hyper_rectangle (triangulation,meshDimensions, min, max);
std::cout << " Number of cells: "
<< triangulation.n_cells()
<< std::endl;
typename Triangulation<dim>::cell_iterator
cell = triangulation.begin (),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int face_number=0;face_number<GeometryInfo<dim>::faces_per_cell;++face_number)
{
if(cell->face(face_number)->at_boundary())
{
if(cell->face(face_number)->center()(0)>0.0)
cell->face(face_number)->set_boundary_id (1);
}
}
}
}
template<int dim>
void Eigen<dim>::deg_of_freedom()
{
dof_handler.distribute_dofs (fe);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
std::cout << "Number of degrees of freedom :"
<< dof_handler.n_dofs()
<<"\ndof of each cell: "
<< dofs_per_cell
<< std::endl;
DoFTools::make_zero_boundary_constraints (dof_handler,0,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());
//Enter the global node x-coordinates into the vector "nodeLocation"
//MappingQ1<dim,dim> mapping;
// std::vector< Point<dim,double> > dof_coords(dof_handler.n_dofs());
// nodeLocation.resize(dof_handler.n_dofs());
// DoFTools::map_dofs_to_support_points<dim,dim>(mapping,dof_handler,dof_coords);
//for(unsigned int i=0; i<dof_coords.size(); i++){
// nodeLocation[i] = dof_coords[i][0];
// }
IndexSet eigenfunction_index_set = dof_handler.locally_owned_dofs ();
eigenfunctions.resize (5);
for (unsigned int i=0; i<eigenfunctions.size (); ++i)
eigenfunctions[i].reinit (eigenfunction_index_set, MPI_COMM_WORLD);
eigenvalues.resize (eigenfunctions.size ());
}
/*
template<int dim>
void Eigen<dim>::define_boundary_conds()
{
//Enter the global node x-coordinates into the vector "nodeLocation"
MappingQ1<dim,dim> mapping;
std::vector< Point<dim,double> > dof_coords(dof_handler.n_dofs());
nodeLocation.resize(dof_handler.n_dofs());
DoFTools::map_dofs_to_support_points<dim,dim>(mapping,dof_handler,dof_coords);
for(unsigned int i=0; i<dof_coords.size(); i++){
nodeLocation[i] = dof_coords[i][0];
//std::cout<<"nodeLocation["<<i<<"]="<< nodeLocation[i]<<std::endl;
}
//Specify boundary condtions (call the function)
// define_boundary_conds();
//const unsigned int totalNodes = dof_handler.n_dofs()/2;
const unsigned int TotalDof = dof_handler.n_dofs();
for(unsigned int globalNode=0; globalNode<TotalDof; globalNode++){
if(nodeLocation[globalNode] == 0){
boundary_values[globalNode] = 0;
}
if(nodeLocation[globalNode] == L){
boundary_values[globalNode] =0;
}
// std::cout<<"boundary_values["<<globalNode<<"]="<<boundary_values[globalNode]<<std::endl;
}
}
*/
template<int dim>
void Eigen<dim>::assembly()
{
const unsigned int no_of_element=triangulation.n_cells();
const unsigned int dofs_per_cell = fe.dofs_per_cell;
FullMatrix<double> globledof(no_of_element,dofs_per_cell);
for(int i=0;i<no_of_element;i++){
for(int j=0;j<dofs_per_cell;j++){
globledof(i,j)=2*i+j;
std::cout << "globledof["<<i<<"][ "<<j<<"]"
<< globledof(i,j)
<< std::endl;
}
}
FullMatrix<double> cell_stiffness_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
//double h_e ;
//std::vector<unsigned int> local_dof_indices (dofs_per_cell);
//std::vector<double> localdof_f_each(dofs_per_cell,0);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
//h_e = nodeLocation[local_dof_indices[1]] - nodeLocation[local_dof_indices[0]];
/*for(int p=0;p<no_of_element;p++){
for (int q=0;q<dofs_per_cell;q++)
{
local_dof_indices.at(q) = globledof(p,q);
std::cout << "localdof_f_each["<<q<<"]"
<< localdof_f_each[q]
<< std::endl;
}*/
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active (),
endc = dof_handler.end ();
for (; cell!=endc; ++cell)
{
//std::cout << "The vector elements are: ";
// for (int i = 0; i < local_dof_indices.size(); i++)
// std::cout <<"globaldof0felement"<< local_dof_indices[i] << "\n ";
cell_stiffness_matrix = 0;
cell_mass_matrix = 0;
double E=2.1e11, a=2, ro=7800, A=30 ,bre=50/1000,hei=5/1000,RAa;
double I=bre*pow(hei,3)/12;
double z=E*I/(2*pow(a,3));
cell_stiffness_matrix (0,0) =z*3;
cell_stiffness_matrix (0,1)=z*3*a;
cell_stiffness_matrix (0,2)=-3*z;
cell_stiffness_matrix (0,3)=z*3*a;
cell_stiffness_matrix (1,0)=3*a*z;
cell_stiffness_matrix (1,1)=4*pow(a,2)*z;
cell_stiffness_matrix (1,2)=-3*a*z;
cell_stiffness_matrix (1,3)=2*pow(a,2)*z;
cell_stiffness_matrix (2,0)=-3*z;
cell_stiffness_matrix (2,1)=-z*3*a;
cell_stiffness_matrix (2,2)=3*z;
cell_stiffness_matrix (2,3)=-3*a*z;
cell_stiffness_matrix (3,0)=3*a*z;
cell_stiffness_matrix (3,1)=2*pow(a,2)*z;
cell_stiffness_matrix (3,2)=-3*a*z;
cell_stiffness_matrix (3,3)=4*pow(a,2)*z;
RAa=ro*A*a/105;
cell_mass_matrix (0,0)=78*RAa;
cell_mass_matrix (0,1)=22*a*RAa;
cell_mass_matrix (0,2)=27*RAa;
cell_mass_matrix (0,3)=-13*a*RAa;
cell_mass_matrix (1,0)=22*a*RAa;
cell_mass_matrix (1,1)=8*pow(a,2)*RAa;
cell_mass_matrix (1,2)=13*a*RAa;
cell_mass_matrix (1,3)=-6*pow(a,2)*RAa;
cell_mass_matrix (2,0)=27*RAa;
cell_mass_matrix (2,1)=13*a*RAa;
cell_mass_matrix (2,2)=78*RAa;
cell_mass_matrix (2,3)=22*a*RAa;
cell_mass_matrix (3,0)=-13*a*RAa;
cell_mass_matrix (3,1)=-6*pow(a,2)*RAa;
cell_mass_matrix (3,2)=22*a*RAa;
cell_mass_matrix (3,3)=8*pow(a,2)*RAa;
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);
}
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>
unsigned int Eigen<dim>::solve ()
{
// We start here, as we normally do, by assigning convergence control we
// want:
SolverControl solver_control (dof_handler.n_dofs(), 1e-9);
SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control);
// Before we actually solve for the eigenfunctions and -values, we have to
// also select which set of eigenvalues to solve for. Lets select those
// eigenvalues and corresponding eigenfunctions with the smallest real
// part (in fact, the problem we solve here is symmetric and so the
// eigenvalues are purely real). After that, we can actually let SLEPc do
// its work:
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 ();
// Finally return the number of iterations it took to converge:
return solver_control.last_step ();
}
template <int dim>
void Eigen<dim>::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>
void Eigen<dim>::run ()
{
generate_mesh(10);
deg_of_freedom();
assembly();
const unsigned int n_iterations = solve ();
std::cout << " Solver converged in " << n_iterations
<< " iterations." << std::endl;
output_results ();
std::cout << std::endl;
for (unsigned int i=0; i<eigenvalues.size(); ++i)
std::cout << " Eigenvalue " << i
<< " : " << eigenvalues[i]
<< std::endl;
// define_boundary_conds();
}
int main (int argc, char **argv)
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
Eigen <1> problem ;
problem.run ();
return 0;
}
