Hello, I think I'm having a similar problem to the OP. I searched a lot, but that's the closest that I could find. I'm trying to work with some meshes from gmsh, so I used the gridin functions to read them.
Triangulation<dim> triangulation;
GridIn<dim> gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("plate4");
gridin.read_msh(f);
The mesh was generated with the physical lines and physical surfaces as
mentioned in step 49.
I even tested if it was actually reading by printing the number of active
elements, and it is.
However, when I'm going to distribute the dofs, it gives this message
An error occurred in line <999> of file
</home/erik/Downloads/dealii-9.0.1/source/dofs/dof_handler.cc> in function
void dealii::DoFHandler<dim, spacedim>::distribute_dofs(const
dealii::FiniteElement<dim, spacedim>&) [with int dim = 2; int spacedim = 2]
The violated condition was:
tria->n_levels() > 0
Additional information:
The Triangulation you are using is empty!
There's some problem that for some reason the information is lost in the
next function when I call
dof_handler.distribute_dofs (fe);
I double checked, and the DoF_handler is already associated with the
triangulation
fe (FE_Q<dim>(order), dim),
dof_handler (triangulation),
It is specially confusing to me, because if I generate a grid with the grid
generator, it works perfectly fine.
Respectfully
E.
Em sábado, 31 de agosto de 2013 07:16:10 UTC-3, Wolfgang Bangerth escreveu:
>
>
> > void CellProb::make_grid ()
> > {
> > Triangulation<2> tria1;
> > GridGenerator::hyper_cube_with_cylindrical_hole (tria1, 0.25, 1.0);
> >
> > Triangulation<2> tria3;
> > GridGenerator::hyper_cube_with_cylindrical_hole (tria3, 0.25, 1.0);
> > GridTools::shift (Point<2>(0,-2), tria3);
> > Triangulation<2> triangulation2;
> > GridGenerator::merge_triangulations (tria1, tria3, triangulation2);
> >
> > mesh_info(triangulation2, "grid-2.eps");
> > std::cout << "Number of active cells: "
> > << triangulation.n_active_cells()
> > << std::endl;
> > std::cout << "Total number of cells: "
> > << triangulation.n_cells()
> > << std::endl;
> >
> > }
> > but I get the following report:
> > ============================ Remaking Makefile.dep
> > ==============debug========= cell.cc -> cell.g.o
> > ============================ Linking cell
> > ============================ Running cell
> > terminate called after throwing an instance of
> > 'dealii::internal::Triangulation::ExcGridHasInvalidCell'
> > what(): --------------------------------------------------------
> > An error occurred in line <1672> of file
> > </usr/local_rwth/sw/deal.II/deal.II-7.3.0/source/grid/tria.cc> in
> function
> > static void
> >
> dealii::internal::Triangulation::Implementation::create_triangulation(const
> > std::vector<dealii::Point<spacedim, double>,
> > std::allocator<dealii::Point<spacedim, double>>> &, const
> > std::vector<dealii::CellData<2>, std::allocator<dealii::CellData<2>>> &,
> const
> > dealii::SubCellData &, dealii::Triangulation<2, spacedim> &) [with
> spacedim = 2]
> > The violated condition was:
> > needed_lines.find(std::make_pair(line_vertices.second,
> > line_vertices.first)) == needed_lines.end()
> > The name and call sequence of the exception was:
> > ExcGridHasInvalidCell(cell)
> > Additional Information:
> > Something went wrong when making cell 9. Read the docs and the source
> code for
> > more information.
> > --------------------------------------------------------
>
> I can reproduce this. Can you please open a bug report with the issue
> tracker,
> including your small testcase that shows the problem?
>
>
> > My second question:
> > I tried grid_2 (from tutorial 49) with the simple version of Poisson's
> > equation and I get the following:
> > ============================ Remaking Makefile.dep
> > ==============debug========= cell.cc -> cell.g.o
> > ============================ Linking cell
> > ============================ Running cell
> > Mesh info:
> > dimension: 2
> > no. of cells: 224
> > boundary indicators: 0(88 times)
> > written to grid-2.eps
> > Number of active cells: 224
> > Total number of cells: 294
> > --------------------------------------------------------
> > An error occurred in line <230> of file
> >
> </usr/local_rwth/sw/deal.II/deal.II-7.3.0/source/dofs/dof_handler_policy.cc>
>
> > in function
> > static unsigned int
> >
> dealii::internal::DoFHandler::Policy::Implementation::distribute_dofs(unsigned
>
>
> > int, unsigned int, dealii::DoFHandler<dim, spacedim> &) [with dim = 2,
> > spacedim = 2]
> > The violated condition was:
> > tria.n_levels() > 0
> > The name and call sequence of the exception was:
> > ExcMessage("Empty triangulation")
> > Additional Information:
> > Empty triangulation
> > Stacktrace:
> > -----------
> > #0 /usr/local_rwth/sw/deal.II/deal.II-7.3.0/lib/libdeal_II.g.so.7.3.0:
> > unsigned int
> > dealii::internal::DoFHandler::Policy::Implementation::distribute_dofs<2,
> > 2>(unsigned int, unsigned int, dealii::DoFHandler<2, 2>&)
> > #1 /usr/local_rwth/sw/deal.II/deal.II-7.3.0/lib/libdeal_II.g.so.7.3.0:
> > dealii::internal::DoFHandler::Policy::Sequential<2,
> > 2>::distribute_dofs(dealii::DoFHandler<2, 2>&) const
> > #2 /usr/local_rwth/sw/deal.II/deal.II-7.3.0/lib/libdeal_II.g.so.7.3.0:
> > dealii::DoFHandler<2, 2>::distribute_dofs(dealii::FiniteElement<2, 2>
> const&)
> > #3 ./cell: CellProb::setup_system()
> > #4 ./cell: CellProb::run()
> > #5 ./cell: main
> > --------------------------------------------------------
> > make: *** [run] Aborted (core dumped)
> > Why is the triangulation empty?
> > I called "triangulation.refine_global (2)" because the documentation
> says
> > "since it will most likely be repopulated soon by the next refinement
> > process", but it did not help.
>
> Can you send us the code you used? This seems odd. Are you sure you have
> associated the dof_handler object with the same triangulation for which
> you
> output the mesh info?
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@math.tamu.edu
> <javascript:>
> www: http://www.math.tamu.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 dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/c55e8e0d-02f0-4647-b767-c37d3903a343%40googlegroups.com.
/*This is a template file for use with 3D finite elements (vector field).
The portions of the code you need to fill in are marked with the comment "//EDIT".
Do not change the name of any existing functions, but feel free
to create additional functions, variables, and constants.
It uses the deal.II FEM library.*/
//Include files
//Data structures and solvers
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/tensor_function.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/sparse_direct.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/base/symmetric_tensor.h>
//Mesh related classes
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.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_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/grid/grid_in.h>
#include <deal.II/grid/manifold_lib.h>
#include <deal.II/grid/grid_out.h>
//Finite element implementation classes
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/precondition.h>
//Standard C++ libraries
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <math.h>
using namespace dealii;
//Define the order of the basis functions (Lagrange polynomials)
//and the order of the quadrature rule globally
const unsigned int order = 1;
const unsigned int quadRule = 2;
template <int dim>
class FEM
{
public:
//Class functions
FEM (); // Class constructor
~FEM(); //Class destructor
//Function to calculate components of the elasticity tensor
double C(unsigned int i,unsigned int j,unsigned int k,unsigned int l);
//Solution steps
void generate_mesh(std::vector<unsigned int> numberOfElements);
void define_boundary_conds();
void setup_system();
void assemble_system();
void solve();
void compute_stress();
void output_results();
//Class objects
Triangulation<dim> triangulation; //mesh
FESystem<dim> fe; //FE element
DoFHandler<dim> dof_handler; //Connectivity matrices
//NEW - deal.II quadrature
QGauss<dim> quadrature_formula; //Quadrature
QGauss<dim-1> face_quadrature_formula; //Face Quadrature
//Data structures
SparsityPattern sparsity_pattern; //Sparse matrix pattern
SparseMatrix<double> K; //Global stiffness matrix - Sparse matrix - used in the solver
Vector<double> D, F; //Global vectors - Solution vector (D) and Global force vector (F)
Table<2,double> dofLocation; //Table of the coordinates of dofs by global dof number
std::map<unsigned int,double> boundary_values; //Map of dirichlet boundary conditions
//solution name array
std::vector<std::string> nodal_solution_names;
std::vector<DataComponentInterpretation::DataComponentInterpretation> nodal_data_component_interpretation;
};
// Class constructor for a vector field
template <int dim>
FEM<dim>::FEM ()
:
fe (FE_Q<dim>(order), dim),
dof_handler (triangulation),
quadrature_formula(quadRule),
face_quadrature_formula(quadRule)
{
//Nodal Solution names - this is for writing the output file
for (unsigned int i=0; i<dim; ++i){
nodal_solution_names.push_back("u");
nodal_data_component_interpretation.push_back(DataComponentInterpretation::component_is_part_of_vector);
}
}
//Class destructor
template <int dim>
FEM<dim>::~FEM (){dof_handler.clear ();}
//Function to calculate the components of the 4th order elasticity tensor
template <int dim>
double FEM<dim>::C(unsigned int i,unsigned int j,unsigned int k,unsigned int l){
double C11,C12,C13,C22,C23,C33,C44,C55,C66,lambda,mu;
C11=8.75*pow(10,9);
C12=1.5*pow(10,8);
C13=0;
C22=8.75*pow(10,9);
C23=0;
C33=0;
C44=0;
C55=0;
C66=4.2*pow(10,9);
lambda=3;
mu=4.3333;
if (i==0 && j==0 && k==0 && l==0)
{
return (C11);
}
else if ((i==0 && j==0 && k==1 && l==1) || (i==1 && j==1 && k==0 && l==0))
{
return (C12);
}
else if ((i==0 && j==0 && k==2 && l==2) || (i==2 && j==2 && k==0 && l==0))
{
return (C13);
}
else if ((i==1 && j==1 && k==1 && l==1))
{
return (C22);
}
else if ((i==2 && j==2 && k==1 && l==1) || (i==1 && j==1 && k==2 && l==2))
{
return (C23);
}
else if ((i==2 && j==2 && k==2 && l==2))
{
return (C33);
}
else if ((i==1 && j==2 && k==1 && l==2) || (i==1 && j==2 && k==2 && l==1) || (i==2 && j==1 && k==1 && l==2) || (i==2 && j==1 && k==2 && l==1))
{
return (C44);
}
else if ((i==0 && j==2 && k==0 && l==2) || (i==0 && j==2 && k==2 && l==0) || (i==2 && j==0 && k==0 && l==2) || (i==2 && j==0 && k==2 && l==0))
{
return (C55);
}
else if ((i==0 && j==1 && k==0 && l==1) || (i==0 && j==1 && k==1 && l==0) || (i==1 && j==0 && k==0 && l==1) || (i==1 && j==0 && k==1 && l==0))
{
return (C66);
}
else
{
return 0;
}
}
//Define the problem domain and generate the mesh
template <int dim>
void FEM<dim>::generate_mesh(std::vector<unsigned int> numberOfElements){
{
Triangulation<dim> triangulation;
GridIn<dim> gridin;
gridin.attach_triangulation(triangulation);
std::ifstream f("plate4");
gridin.read_msh(f);
triangulation.refine_global(1);
std::cout << " Number of active elems: " << triangulation.n_active_cells() << std::endl;
}
}
//Specify the Dirichlet boundary conditions
template <int dim>
void FEM<dim>::define_boundary_conds(){
//EDIT - Define the Dirichlet boundary conditions.
/*Note: this will be very similiar to the define_boundary_conds function
in the HW2 template. You will loop over all degrees of freedom and use "dofLocation"
to check if the dof is on the boundary with a Dirichlet condition.
(Aside: Since we now have more than 1 dof per node, it is possible to apply a different Dirichlet
condition on each of the 3 dofs on the same node. Note that this is NOT the case for the current
assignment. But if you wanted to do it, you would need to check
the nodal dof (0, 1, or 2) as well as the location. For example, if you wanted to fix displacements
only in the x-direction at x=0, you would have an condition such as this:
(dofLocation[globalDOF][0] == 0 && nodalDOF == 0)
Note that you can get the nodal dof from the global dof using the "modulo" operator,
i.e nodalDOF = globalDOF%dim. Here "%" gives you the remainder when dividing globalDOF by dim.)
Add the global dof number and the specified value (displacement in this problem)
to the boundary values map, something like this:
boundary_values[globalDOFIndex] = dirichletDisplacementValue
Note that "dofLocation" is now a Table of degree of freedom locations, not just node locations,
so, for example, dofs 0, 1, and 2 correspond to the same node, so that have the same coordinates.
The row index is the global dof number; the column index refers to the x, y, or z component (0, 1, or 2 for 3D).
e.g. dofLocation[7][2] is the z-coordinate of global dof 7*/
const unsigned int totalDOFs = dof_handler.n_dofs(); //Total number of degrees of freedom
for(unsigned int globalDOF=0; globalDOF<totalDOFs; globalDOF++){
unsigned int nodalDOF = globalDOF%dim;
if(dofLocation[globalDOF][0]+1 == 0){
boundary_values[globalDOF] = 0;
}
if(dofLocation[globalDOF][1]+1 == 0 ){
boundary_values[globalDOF] = 0;
}
// if(dofLocation[globalDOF][0] == 1 && nodalDOF==0){ //Aplica deslocamento prescrito na direção x na face x = 1
// boundary_values[globalDOF] = 1;
// }
}
}
//Setup data structures (sparse matrix, vectors)
template <int dim>
void FEM<dim>::setup_system(){
std::cout << " Number of elems: " << triangulation.n_cells() << std::endl;
std::cout << " Number of active elems: " << triangulation.n_active_cells() << std::endl;
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl;
//Let deal.II organize degrees of freedom
dof_handler.distribute_dofs (fe);
//Get a vector of global degree-of-freedom x-coordinates
MappingQ1<dim,dim> mapping;
std::vector< Point<dim,double> > dof_coords(dof_handler.n_dofs());
dofLocation.reinit(dof_handler.n_dofs(),dim);
DoFTools::map_dofs_to_support_points<dim,dim>(mapping,dof_handler,dof_coords);
for(unsigned int i=0; i<dof_coords.size(); i++){
for(unsigned int j=0; j<dim; j++){
dofLocation[i][j] = dof_coords[i][j];
}
}
//Specify boundary condtions (call the function)
define_boundary_conds();
//Define the size of the global matrices and vectors
sparsity_pattern.reinit (dof_handler.n_dofs(),
dof_handler.n_dofs(),
dof_handler.max_couplings_between_dofs());
DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
sparsity_pattern.compress();
K.reinit (sparsity_pattern);
D.reinit(dof_handler.n_dofs());
F.reinit(dof_handler.n_dofs());
//Just some notes...
std::cout << " Number of active elems: " << triangulation.n_active_cells() << std::endl;
std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl;
}
//Form elmental vectors and matrices and assemble to the global vector (F) and matrix (K)
template <int dim>
void FEM<dim>::assemble_system(){
/*NEW - deal.II basis functions, etc. The third input values (after fe and quadrature_formula)
specify what information we want to be updated. For fe_values, we need the basis function values,
basis function gradients, and det(Jacobian) times the quadrature weights (JxW). For fe_face_values,
we need the basis function values, the value of x at the quadrature points, and JxW.*/
//For volume integration/quadrature points
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values |
update_gradients |
update_quadrature_points |
update_JxW_values);
//For surface integration/quadrature points
FEFaceValues<dim> fe_face_values (fe,
face_quadrature_formula,
update_values |
update_quadrature_points |
update_JxW_values);
K=0; F=0;
const unsigned int dofs_per_elem = fe.dofs_per_cell; //This gives you dofs per element
const unsigned int nodes_per_elem = fe.dofs_per_cell/dim;
const unsigned int num_quad_pts = quadrature_formula.size(); //Total number of quad points in the element
const unsigned int num_face_quad_pts = face_quadrature_formula.size(); //Total number of quad points in the face
const unsigned int faces_per_elem = GeometryInfo<dim>::faces_per_cell;
FullMatrix<double> Klocal (dofs_per_elem, dofs_per_elem);
Vector<double> Flocal (dofs_per_elem);
std::vector<unsigned int> local_dof_indices (dofs_per_elem); //This relates local dof numbering to global dof numbering
//loop over elements
typename DoFHandler<dim>::active_cell_iterator elem = dof_handler.begin_active(),
endc = dof_handler.end();
for (;elem!=endc; ++elem){
//Update fe_values for the current element
fe_values.reinit(elem);
//Retrieve the effective "connectivity matrix" for this element
elem->get_dof_indices (local_dof_indices);
/*Global Assembly -
Get the current Flocal and Klocal from the functions you wrote above
and populate F_assembly and K_assembly using local_dof_indices to relate local and global DOFs.*/
Klocal = 0.;
//Loop over local DOFs and quadrature points to populate Klocal
//Note that all quadrature points are included in this single loop
for (unsigned int q=0; q<num_quad_pts; ++q){
// double z = fe_values.quadrature_point(q)[2]; //y-coordinate at the current surface quad. point
// std::cout << "z " << z << std::endl;
//evaluate elemental stiffness matrix, K^{AB}_{ik} = \integral N^A_{,j}*C_{ijkl}*N^B_{,l} dV
for (unsigned int A=0; A<nodes_per_elem; A++) { //Loop over nodes
for(unsigned int i=0; i<dim; i++){ //Loop over nodal dofs
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++){
/*//EDIT - You need to define Klocal here. Note that the indices of Klocal are the element dof numbers (0 through 23),
which you can caluclate from the element node numbers (0 through 8) and the nodal dofs (0 through 2).
You'll need the following information:
basis gradient vector: fe_values.shape_grad(elementDOF,q), where elementDOF is dim*A+i or dim*B+k
NOTE: this is the gradient with respect to the real domain (not the bi-unit domain)
elasticity tensor: use the function C(i,j,k,l)
det(J) times the total quadrature weight: fe_values.JxW(q)*/
Klocal[dim*A+i][dim*B+k]+=fe_values.shape_grad(dim*A+i,q)[j]*C(i,j,k,l)*
fe_values.shape_grad(dim*B+k,q)[l]*fe_values.JxW(q);
}
}
}
}
}
}
}
Flocal = 0.;
//Loop over faces (for Neumann BCs), local DOFs and quadrature points to populate Flocal.
//If you had a forcing function, you would need to use FEValues here as well to integrate over the volume.
//Add Neumann boundary conditions here in Flocal by integrating over the appropriate surface
Vector<double> h(dim); h=0.;
for (unsigned int f=0; f < faces_per_elem; f++){
//Update fe_face_values from current element and face
fe_face_values.reinit (elem, f);
/*elem->face(f)->center() gives a position vector (in the real domain) of the center point on face f
of the current element. We can use it to see if we are at the Neumann boundary, x_3 = 1.*/
if(elem->face(f)->center()[0]+0.5 >= 1){
// if(elem->face(f)->center()[2] == 1){
//To integrate over this face, loop over all face quadrature points with this single loop
for (unsigned int q=0; q<num_face_quad_pts; ++q){
//EDIT - define the value of the traction vector, h
// h[2]=1.e9*x;
h[0]=pow(10,5);
// h[0]=1.;
for (unsigned int A=0; A<nodes_per_elem; A++){ //loop over all element nodes
for(unsigned int i=0; i<dim; i++){ //loop over nodal dofs
/*//EDIT - define Flocal. Again, the indices of Flocal are the element dof numbers (0 through 23).
Evaluate the basis functions using the elementDOF: fe_face_values.shape_value(elementDOF,q)
Note that we are looping over all element dofs, not just those on the Neumann face. However,
the face quadrature points are only on the Neumann face, so we are indeed doing a surface integral.
For det(J) times the total quadrature weight: fe_face_values.JxW(q)*/
Flocal[dim*A+i]+=h[i]*fe_face_values.shape_value(dim*A+i,q)*fe_face_values.JxW(q);
}
}
}
}
if(elem->face(f)->center()[0]+0.5 == 0){
// if(elem->face(f)->center()[2] == 1){
//To integrate over this face, loop over all face quadrature points with this single loop
for (unsigned int q=0; q<num_face_quad_pts; ++q){
//EDIT - define the value of the traction vector, h
// h[2]=1.e9*x;
h[0]=-pow(10,5);
// h[1]=1.;
for (unsigned int A=0; A<nodes_per_elem; A++){ //loop over all element nodes
for(unsigned int i=0; i<dim; i++){ //loop over nodal dofs
/*//EDIT - define Flocal. Again, the indices of Flocal are the element dof numbers (0 through 23).
Evaluate the basis functions using the elementDOF: fe_face_values.shape_value(elementDOF,q)
Note that we are looping over all element dofs, not just those on the Neumann face. However,
the face quadrature points are only on the Neumann face, so we are indeed doing a surface integral.
For det(J) times the total quadrature weight: fe_face_values.JxW(q)*/
Flocal[dim*A+i]+=h[i]*fe_face_values.shape_value(dim*A+i,q)*fe_face_values.JxW(q);
}
}
}
}
if(elem->face(f)->center()[1]+0.5 == 0){
// if(elem->face(f)->center()[2] == 1){
//To integrate over this face, loop over all face quadrature points with this single loop
for (unsigned int q=0; q<num_face_quad_pts; ++q){
//EDIT - define the value of the traction vector, h
// h[2]=1.e9*x;
h[1]=0.9*pow(10,5);
// h[1]=1.;
for (unsigned int A=0; A<nodes_per_elem; A++){ //loop over all element nodes
for(unsigned int i=0; i<dim; i++){ //loop over nodal dofs
/*//EDIT - define Flocal. Again, the indices of Flocal are the element dof numbers (0 through 23).
Evaluate the basis functions using the elementDOF: fe_face_values.shape_value(elementDOF,q)
Note that we are looping over all element dofs, not just those on the Neumann face. However,
the face quadrature points are only on the Neumann face, so we are indeed doing a surface integral.
For det(J) times the total quadrature weight: fe_face_values.JxW(q)*/
Flocal[dim*A+i]+=h[i]*fe_face_values.shape_value(dim*A+i,q)*fe_face_values.JxW(q);
}
}
}
}
if(elem->face(f)->center()[1]+0.5 == 1){
// if(elem->face(f)->center()[2] == 1){
//To integrate over this face, loop over all face quadrature points with this single loop
for (unsigned int q=0; q<num_face_quad_pts; ++q){
//EDIT - define the value of the traction vector, h
// h[2]=1.e9*x;
h[1]=-0.9*pow(10,5);
// h[1]=1.;
for (unsigned int A=0; A<nodes_per_elem; A++){ //loop over all element nodes
for(unsigned int i=0; i<dim; i++){ //loop over nodal dofs
/*//EDIT - define Flocal. Again, the indices of Flocal are the element dof numbers (0 through 23).
Evaluate the basis functions using the elementDOF: fe_face_values.shape_value(elementDOF,q)
Note that we are looping over all element dofs, not just those on the Neumann face. However,
the face quadrature points are only on the Neumann face, so we are indeed doing a surface integral.
For det(J) times the total quadrature weight: fe_face_values.JxW(q)*/
Flocal[dim*A+i]+=h[i]*fe_face_values.shape_value(dim*A+i,q)*fe_face_values.JxW(q);
}
}
}
}
}
//Assemble local K and F into global K and F
for(unsigned int i=0; i<dofs_per_elem; i++){
//EDIT - Assemble F from Flocal (you can look at HW2)
F[local_dof_indices[i]]+=Flocal[i];
for(unsigned int j=0; j<dofs_per_elem; j++){
//EDIT - Assemble K from Klocal (you can look at HW2)
K.add(local_dof_indices[i],local_dof_indices[j],Klocal[i][j]);
}
}
}
//Let deal.II apply Dirichlet conditions WITHOUT modifying the size of K and F global
MatrixTools::apply_boundary_values (boundary_values, K, D, F, false);
}
//Solve for D in KD=F
template <int dim>
void FEM<dim>::solve(){
//Solve for D
SolverControl solver_control(1000000, 1e-12);
SolverCG<> solver(solver_control);
solver.solve(K, D, F, PreconditionIdentity());
// std::cout << "D " << D << std::endl;
}
template <int dim>
void FEM<dim>::compute_stress(){
std::ofstream strain;
strain.open ("strain.txt");
// QTrapez<dim> vertex_quadrature;
// std::vector<Tensor<1, dim>> solution_gradients(vertex_quadrature.size());
FEValues<dim> fe_values(fe,
quadrature_formula,
update_values |
update_gradients |
update_quadrature_points |
update_JxW_values);
// fe_values.get_function_gradients(D, solution_gradients);
const unsigned int dofs_per_elem = fe.dofs_per_cell; //This gives you dofs per element
std::vector<unsigned int> local_dof_indices (dofs_per_elem);
const unsigned int num_quad_pts = quadrature_formula.size(); //Total number of quad points in the element
double x, y, strain1, strain2, strain3,E1,E2,G12,null12,null21;
Vector<double> stress_vector(3);
Vector<double> strain_vector(3);
E1=14*pow(10,9),E2=3.5*pow(10,9),G12=4.2*pow(10,9),null21=null12*E2/E1;
SymmetricTensor<2, 3> strain_stress;
strain_stress[0][0] = E1/(1-null12*null21);
strain_stress[1][1] = E2/(1-null12*null21);
strain_stress[2][2] = G12;
strain_stress[0][1] = null21*E1/(1-null21*null12);
strain_stress[0][2] = 0;
strain_stress[1][2] = 0;
// std::cout << "strain_stress " << strain_stress << std::endl;
//loop over elements
typename DoFHandler<dim>::active_cell_iterator elem = dof_handler.begin_active (),
endc = dof_handler.end();
for (;elem!=endc; ++elem){
fe_values.reinit(elem);
//Retrieve the effective "connectivity matrix" for this element
elem->get_dof_indices (local_dof_indices);
for(unsigned int q=0; q<num_quad_pts; q++){
strain1=0;
strain2=0;
strain3=0;
x = fe_values.quadrature_point(q)[0]+0.5; //x-coordinate at the current surface quad. point
y = fe_values.quadrature_point(q)[1]+0.5; //y-coordinate at the current surface quad. point
// std::cout << "x " << x << std::endl;
// std::cout << "y " << y << std::endl;
SymmetricTensor<2, dim> tmp;
//Find the values of x and u_h (the finite element solution) at the quadrature points
for(unsigned int B=0; B<dofs_per_elem; B++){
for (unsigned int i = 0; i < dim; ++i)
tmp[i][i] = fe_values.shape_grad_component(B, q, i)[i];
for (unsigned int i = 0; i < dim; ++i)
for (unsigned int j = i + 1; j < dim; ++j)
tmp[i][j] =
(fe_values.shape_grad_component(B, q, i)[j] +
fe_values.shape_grad_component(B, q, j)[i]) /
2;
if (B%2==0){
strain1 += D[local_dof_indices[B]]*tmp[0][0];
strain3 += D[local_dof_indices[B]]*tmp[0][1];
}
else{
strain2 += D[local_dof_indices[B]]*tmp[1][1];
strain3 += D[local_dof_indices[B]]*tmp[0][1];
}}
strain << x << " ";
strain << y << " ";
strain << strain1 << " ";
strain << strain2 << " ";
strain << strain3 << std::endl;
strain_vector[0]=strain1;
strain_vector[1]=strain2;
strain_vector[2]=strain3;
for (unsigned int i = 0; i < 3; ++i){
for (unsigned int j = 0; j < 3; ++j){
stress_vector[i]+=strain_stress[i][j]*strain_vector[j];
}
}
// std::cout << "stress_vector " << stress_vector << std::endl;
// std::cout << "strain1 " << strain1 << std::endl;
// std::cout << "strain2 " << strain2 << std::endl;
// std::cout << "strain3 " << strain3 << std::endl;
}
}
strain.close();
}
//Output results
template <int dim>
void FEM<dim>::output_results (){
//Write results to VTK file
std::ofstream output1 ("solution.vtk");
DataOut<dim> data_out; data_out.attach_dof_handler (dof_handler);
//Add nodal DOF data
data_out.add_data_vector (D,
nodal_solution_names,
DataOut<dim>::type_dof_data,
nodal_data_component_interpretation);
data_out.build_patches();
data_out.write_vtk(output1);
output1.close();
}
/*This is a skeleton code file for use with the Finite Element Method for Problems in Physics.
It uses the deal.II FEM library, dealii.org*/
//Include files
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include "FEM3.h"
//#include "writeSolutions.h"
using namespace dealii;
//The main program, using the FEM class
int main (){
try{
deallog.depth_console (0);
const int dimension = 2;
//NOTE: This is where you define the number of elements in the mesh
std::vector<unsigned int> num_of_elems(dimension);
num_of_elems[0] = 5;
num_of_elems[1] = 5;
// num_of_elems[2] = 10; //For example, a 10 x 10 x 10 element mesh
FEM<dimension> problemObject;
problemObject.generate_mesh(num_of_elems);
problemObject.setup_system();
problemObject.assemble_system();
problemObject.solve();
problemObject.compute_stress();
problemObject.output_results();
//write solutions to h5 file
// char tag[21];
// sprintf(tag, "CA3");
// writeSolutionsToFileCA3(problemObject.D, tag);
}
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;
}
plate4
Description: Binary data
plate4.geo
Description: Binary data
