Hi all,
Since the sequential code is solved by SolverCG+SSOR and the parallel one
isn't but by Bicgstab and GMRES, It seems to me that I have ruined the
symmetry of my system_matrix through parallelization.
The assembly and make_constraint parts which contribute to the creation and
features of the system_matrix in sequential and parallel code are identical
except in a couple of added lines in the parallel which are highlighted
bellow:
///////////////////////////////////////////////////////// make_constraints
template <int dim>
void Solid<dim>::make_constraints(const int &it_nr)
{
if (it_nr > 1)
return;
constraints.clear();
constraints.reinit (locally_relevant_dofs);
const bool apply_dirichlet_bc = (it_nr == 0);
const FEValuesExtractors::Scalar x_displacement(0);
const FEValuesExtractors::Scalar y_displacement(1);
const FEValuesExtractors::Scalar z_displacement(2);
{
const int boundary_id = 0;
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
{
const int boundary_id = 1;
const int timestep = time.get_timestep();
if (apply_dirichlet_bc == true)
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
BoundaryDisplacement<dim>(0, timestep),
constraints,
fe.component_mask(x_displacement));
else
VectorTools::interpolate_boundary_values(dof_handler,
boundary_id,
ZeroFunction<dim>(dim),
constraints,
fe.component_mask(x_displacement));
}
constraints.close();
}
////////////////////////////////////////////////////////////////////
assembly
template <int dim>
void Solid<dim>::assemble_system ()
{
tangent_matrix = 0;
system_rhs = 0;
FEValues<dim> fe_values (fe, qf_cell,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FEFaceValues<dim> fe_face_values (fe, qf_face,
update_values |
update_quadrature_points |
update_normal_vectors |
update_JxW_values);
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs (dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
std::vector<double> Nx(dofs_per_cell);
std::vector<Tensor<2, dim> > grad_Nx(dofs_per_cell);
std::vector<SymmetricTensor<2, dim> > symm_grad_Nx(dofs_per_cell);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
{
fe_values.reinit (cell);
cell_matrix = 0;
cell_rhs = 0;
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
{
const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv();
const Tensor<2, dim> tau = lqph[q_point].get_tau();
const SymmetricTensor<2, dim> symm_tau =
lqph[q_point].get_tau();
const SymmetricTensor<4, dim> Jc = lqph[q_point].get_Jc();
const double JxW = fe_values.JxW(q_point);
for (unsigned int k=0; k<dofs_per_cell; ++k)
{
grad_Nx[k] = fe_values[u_fe].gradient(k, q_point) * F_inv;
symm_grad_Nx[k] = symmetrize(grad_Nx[k]);
}
for (unsigned int i=0; i<dofs_per_cell; ++i)
{
const unsigned int component_i =
fe.system_to_component_index(i).first;
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
const unsigned int component_j =
fe.system_to_component_index(j).first;
cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material
contribution:
* symm_grad_Nx[j] * JxW;
if (component_i == component_j) // geometrical stress
contribution
cell_matrix(i, j) += grad_Nx[i][component_i] * tau
* grad_Nx[j][component_j] * JxW;
}
cell_rhs(i) -= symm_grad_Nx[i] * symm_tau * JxW;
}
}
cell->get_dof_indices (local_dof_indices);
constraints.distribute_local_to_global (cell_matrix,
cell_rhs,
local_dof_indices,
tangent_matrix,
system_rhs);
}
tangent_matrix.compress (VectorOperation::add);
system_rhs.compress (VectorOperation::add);
}
///////////////////////////////////////////////////////////////
I was wondering how and where I may have destroyed the symmetry condition.
Is there any other probable function that contributes to the symmetry of my
system_matrix.
In addition I don't know how to plot the system_matrix components to
compare system_matrix of sequential and parallel codes.
Thanks,
--
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.
