Hi,
Sorry to disturb again, I have a problem. :)
The following line gives me a headache:
std::cout << cell->face(face_number)->boundary_id() << std::endl;
... since I cannot see anything. It outputs strange symbols such as � or
some 0001 cube.
My aim is to somehow write my RHS for distributed loads or point loads
(Neumann BCs) as well as prescribed displacements (Dirichlet BCs).
In order to do so, I took step-7 and step-17 as a guide, because it should
work with MPI also.
My RHS class definition looks like this:
template<int dim>
class RightHandSide: public Function<dim>
{
public:
RightHandSide () :
Function<dim>()
{
}
virtual double value (const Point<dim> &p, const unsigned int c = 0) const;
};
template<int dim>
double RightHandSide<dim>::value (const Point<dim> &p, const unsigned int)
const
{
if ( std::fabs(p(1)) + 1.0 < 1.e-11 ) // y = -5.0 is one of the domain
boundaries return 1.0;
return 1.;
else
return 0.;
}
The value function is just dummy, to see if I even get any values inside my
RHS vector.
My assembly function:
template<int dim>
void SolidMechanics<dim>::assemble_system ()
{
// Timer for "assembly"
TimerOutput::Scope t(computing_timer, "assembly");
// Initialize stiffness matrix and external forces
system_matrix = 0;
system_rhs = 0;
// Initialize quadrature integration
QGauss<dim> quadrature_formula(2);
QGauss<dim - 1> face_quadrature_formula(2);
// Compute FE values (e.g. shape functions | shape function gradients |
quadrature points | Jacobian determinants x quadrature weights)
FEValues<dim> fe_values(fe, quadrature_formula, update_values |
update_gradients | update_quadrature_points | update_JxW_values |
update_jacobians | update_inverse_jacobians);
FEFaceValues<dim> fe_face_values(fe, face_quadrature_formula, update_values
| update_quadrature_points | update_normal_vectors | update_JxW_values);
// Initialization of variables
const unsigned int dofs_per_cell = fe.dofs_per_cell; // Number of degrees
of freedom for each cell
const unsigned int n_q_points = quadrature_formula.size(); // Number of
quadrature points (interior)
const unsigned int n_face_q_points = face_quadrature_formula.size(); //
Number of quadrature points (boundary)
const unsigned int vertices_per_cell =
GeometryInfo<dim>::vertices_per_cell; // Number of vertices for each cell
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell); // Stiffness
matrix
Vector<double> cell_rhs(dofs_per_cell); // External force vector
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell); //
Global number of DoFs
const RightHandSide<dim> right_hand_side;
std::vector<double> rhs_values(n_face_q_points);
// Assemble system
typename DoFHandler<dim>::active_cell_iterator cell =
dof_handler.begin_active(), endc = dof_handler.end();
for (; cell != endc; ++cell) // loop over all cells
if ( cell->is_locally_owned() ) // only loop over locally owned cells (~
current processor)
{
cell_matrix = 0;
cell_rhs = 0;
fe_values.reinit(cell); // re-initialize FE values for current cell
// COMPUTE: Stiffness matrix
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) // loop
over integration/quadrature points (GPs)
{
const SymmetricTensor<2, dim> eps_phi_i = Tensors::get_strain(fe_values, i,
q_point), eps_phi_j = Tensors::get_strain(fe_values, j, q_point);
cell_matrix(i, j) += (eps_phi_i * elasticity_tensor * eps_phi_j *
fe_values.JxW(q_point));
// print_tensor(eps_phi_i, "STRAIN (LEFT)");
// print_tensor(eps_phi_j, "STRAIN (RIGHT)");
}
}
}
// COMPUTE: Internal force vector (NEUMANN boundary)
for (unsigned int face_number = 0; face_number <
GeometryInfo<dim>::faces_per_cell; ++face_number)
{
std::cout << cell->face(face_number)->boundary_id() << std::endl;
if ( cell->face(face_number)->at_boundary() &&
(cell->face(face_number)->boundary_id() == 3) )
{
fe_face_values.reinit(cell, face_number);
std::cout << "Boundary 3 is selected.
============================================" << std::endl;
right_hand_side.value_list(fe_face_values.get_quadrature_points(),
rhs_values);
for (unsigned int q_point = 0; q_point < n_face_q_points; ++q_point)
{
const Tensor<1, dim> neumann_value = (rhs_values[q_point]) *
fe_face_values.normal_vector(q_point);
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
const unsigned int component_i = fe.system_to_component_index(i).first;
if ( component_i < dim )
cell_rhs(i) += (neumann_value[component_i] * fe_face_values.shape_value(i,
q_point) * fe_face_values.JxW(q_point));
}
}
}
}
// Global numbers of degrees of freedom on this cell
cell->get_dof_indices(local_dof_indices);
// Copy the local matrix and vector on each cell into the global system
constraints.distribute_local_to_global(cell_matrix, cell_rhs,
local_dof_indices, system_matrix, system_rhs);
print_vector(cell_rhs, "\nELEMENT INTERNAL FORCE VECTOR\n");
// check_values("internal_forces", fe_values, dofs_per_cell, n_q_points,
vertices_per_cell, elasticity_tensor, cell_matrix, cell_rhs); // check
values
}
According to step-7 I also use this algorithm which I think should set the
boundary ID:
// Geometry: Simple cube
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(1);
typename parallel::distributed::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 ( (std::fabs(cell->face(face_number)->center()(0) - (-1)) < 1e-12) ||
(std::fabs(cell->face(face_number)->center()(1) - (-1)) < 1e-12) )
cell->face(face_number)->set_boundary_id(1);
I already tried to learn from the recent posts, e.g.
https://groups.google.com/forum/#!searchin/dealii/value_list%7Csort:relevance/dealii/zGhX17espQI/wiF41ZydDwAJ
or
https://groups.google.com/forum/#!searchin/dealii/$20$20pressure_boundary_values|sort:relevance/dealii/S2SR6cgncs8/En56YKJACgAJ
..., but without success.
The RHS vector stays zero. No entries. My question is just what would be
the correct "value" or "vector_value_list" function to define in order to
apply traction or pressure at the upper/top boundary of a simple cube. But
prior to that I have to see why my boundary ID is not correctly found.
Your help would be greatly appreciated, thank you in advance.
Kind regards,
S. A. Mohseni
--
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.
