Hi all,
Following step-45, I am going to implement periodic boundary condition on
my code for a Thermo-elastic problem, in which the thermal and elastic
parts are solved separately (I have two different dof_handelers for
temprature and displacement fields without any block structure). I am
wondering how to apply periodic constrains, made in setup_system, in the
thermal field assemble_system for which I don't use
"constraints.distribute_local_to_global ()" function. Thermal
Assemble_system is the following:
template <int dim>
void Solid<dim>::assemble_system_eta ()
{
system_matrix_eta = 0;
system_rhs_eta = 0;
mass_matrix = 0;
laplace_matrix = 0;
nl_matrix = 0;
nl_term = 0;
FEValues<dim> fe_values_eta (fe_eta, qf_cell_eta,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FullMatrix<double> cell_mass_matrix (dofs_per_cell_eta,
dofs_per_cell_eta);
FullMatrix<double> cell_laplace_matrix (dofs_per_cell_eta,
dofs_per_cell_eta);
FullMatrix<double> cell_nl_matrix (dofs_per_cell_eta,
dofs_per_cell_eta);
Vector <double> cell_nl_term (dofs_per_cell_eta);
std::vector<types::global_dof_index>
local_dof_indices(dofs_per_cell_eta);
std::vector<double> phi (dofs_per_cell_eta);
std::vector<Tensor<1,dim> > grad_phi (dofs_per_cell_eta);
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler_eta.begin_active(),
endc = dof_handler_eta.end();
for (; cell!=endc; ++cell)
if (cell->is_locally_owned())
{
fe_values_eta.reinit(cell);
cell_mass_matrix = 0;
cell_laplace_matrix = 0;
cell_nl_matrix = 0;
cell_nl_term = 0;
PointHistory<dim> *lqph =
reinterpret_cast<PointHistory<dim>*>(cell->user_pointer());
for (unsigned int q=0; q<n_q_points_eta; ++q)
{
const double driving_force = lqph[q].get_driving_force();
const double deri_driving_force = lqph[q].get_deri_driving_force();
for (unsigned int k=0; k<dofs_per_cell_eta; ++k)
{
phi[k] = fe_values_eta.shape_value (k, q);
grad_phi[k] = fe_values_eta.shape_grad (k, q);
}
for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
{
for (unsigned int j=0; j<dofs_per_cell_eta; ++j)
{
cell_mass_matrix(i,j) += (phi[i] * phi[j]) *
fe_values_eta.JxW(q);
cell_laplace_matrix(i,j) += (grad_phi[i] * grad_phi[j]) *
fe_values_eta.JxW(q);
cell_nl_matrix(i,j) += deri_driving_force * phi[i] * phi[j]
* fe_values_eta.JxW(q);
}
cell_nl_term(i) += driving_force * phi[i] *
fe_values_eta.JxW (q);
}
}
cell->get_dof_indices (local_dof_indices);
for (unsigned int i=0; i<dofs_per_cell_eta; ++i)
{
for (unsigned int j=0; j<dofs_per_cell_eta; ++j)
{
mass_matrix. add (local_dof_indices[i],
local_dof_indices[j],
cell_mass_matrix(i,j));
laplace_matrix.add (local_dof_indices[i],
local_dof_indices[j],
cell_laplace_matrix(i,j));
nl_matrix. add (local_dof_indices[i],
local_dof_indices[j],
cell_nl_matrix(i,j));
}
nl_term (local_dof_indices[i]) += cell_nl_term(i);
}
}
mass_matrix.compress (VectorOperation::add);
laplace_matrix.compress (VectorOperation::add);
nl_matrix.compress (VectorOperation::add);
nl_term.compress (VectorOperation::add);
TrilinosWrappers::MPI::Vector temp_solution_eta
(locally_owned_dofs_eta, mpi_communicator);
TrilinosWrappers::MPI::Vector temp_old_solution_eta
(locally_owned_dofs_eta, mpi_communicator);
temp_solution_eta = solution_eta;
temp_old_solution_eta = old_solution_eta;
system_matrix_eta.copy_from (mass_matrix);
system_matrix_eta.add ( parameters.L * parameters.beta *
parameters.delta_t * parameters.thetan, laplace_matrix);
system_matrix_eta.add ( -parameters.L * parameters.delta_t *
parameters.thetan, nl_matrix);
tmp_matrix.copy_from (mass_matrix);
tmp_matrix.add ( parameters.L * parameters.beta * parameters.delta_t
* parameters.thetan, laplace_matrix);
tmp_matrix.vmult (tmp_vector, temp_solution_eta);
system_rhs_eta += tmp_vector;
tmp_matrix.add( -parameters.L * parameters.beta * parameters.delta_t,
laplace_matrix);
tmp_matrix.vmult (tmp_vector, temp_old_solution_eta);
system_rhs_eta -= tmp_vector;
system_rhs_eta.add (-parameters.L * parameters.delta_t, nl_term);
system_rhs_eta *= -1;
system_matrix_eta.compress (VectorOperation::add);
system_rhs_eta.compress (VectorOperation::add);
}
I guessed perhaps I can replace the red part with the following,
constraints.distribute_local_to_global (cell_mass_matrix,
local_dof_indices,
mass_matrix);
constraints.distribute_local_to_global (cell_laplace_matrix,
local_dof_indices,
laplace_matrix);
constraints.distribute_local_to_global (cell_nl_matrix,
local_dof_indices,
nl_matrix);
constraints.distribute_local_to_global (cell_nl_term,
local_dof_indices,
nl_term);
but the newton_rophson solver for thermal part diverged when I did so!!
Any hint would be appreciated.
Thanks,
Hamed
--
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.
