I still have some problems trying to understand when to use multiple
DoFHandler-variables/FE_Q-variables, and when to use only one. I am trying
to expand example 52 to multiple coupled equations. My initial approach is
to expand the functions from (here assemble_system)
FEValues<2> fe_values(fe, quadrature_formula,
update_values | update_gradients |
update_JxW_values);
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
FullMatrix<double> cell_mass_matrix (dofs_per_cell, dofs_per_cell);
std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
DoFHandler<2>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
cell_matrix = 0.;
cell_mass_matrix = 0.;
fe_values.reinit (cell);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_matrix(i,j) += ((-diffusion_coefficient *
fe_values.shape_grad(i,q_point) *
fe_values.shape_grad(j,q_point)
- absorption_cross_section *
fe_values.shape_value(i,q_point) *
fe_values.shape_value(j,q_point)) *
fe_values.JxW(q_point));
cell_mass_matrix(i,j) += fe_values.shape_value(i,q_point) *
fe_values.shape_value(j,q_point) *
fe_values.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices);
constraint_matrix.distribute_local_to_global(cell_matrix,
local_dof_indices,system_matrix);
constraint_matrix.distribute_local_to_global(cell_mass_matrix,
local_dof_indices,mass_matrix);
}
to
const QGauss<dim> quadrature_formula(fe_degree+1);
FEValues<dim> fe_values_N(fe, quadrature_formula,
update_values | update_gradients |
update_JxW_values),
fe_values_TE(fe, quadrature_formula, update_values |
update_gradients | update_JxW_values);
const unsigned int dofs_per_cell_N = fe.dofs_per_cell;
const unsigned int dofs_per_cell_TE = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
FullMatrix<double> cell_matrix_N (dofs_per_cell_N, dofs_per_cell_N),
cell_matrix_TE(dofs_per_cell_TE, dofs_per_cell_TE);
FullMatrix<double> cell_mass_matrix_N (dofs_per_cell_N,
dofs_per_cell_N), cell_mass_matrix_TE (dofs_per_cell_TE, dofs_per_cell_TE);
std::vector<types::global_dof_index> local_dof_indices_N (
dofs_per_cell_N), local_dof_indices_TE (dofs_per_cell_TE);
for (auto cell = dof_handler.begin_active(); cell!=dof_handler.end
(); ++cell)
{
cell_matrix_N = 0.;
cell_matrix_TE = 0.;
cell_mass_matrix_N = 0.;
cell_mass_matrix_TE = 0.;
fe_values_N.reinit (cell);
fe_values_TE.reinit (cell);
for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
for (unsigned int i=0; i<dofs_per_cell; ++i)
for (unsigned int j=0; j<dofs_per_cell; ++j)
{
cell_matrix_N(i,j) += ((-diffusion_coefficient_N(
fe_values_N.shape_value(i, q_point)*fe_values_N.shape_value(j, q_point)*
fe_values_N.JxW(q_point)) *
fe_values_N.shape_grad(i,
q_point) *
fe_values_N.shape_grad(j,
q_point)
- absorption_cross_section_N(
fe_values_N.shape_value(i, q_point)*fe_values_N.shape_value(j, q_point)*
fe_values_N.JxW(q_point)) *
fe_values_N.shape_value(i,
q_point) *
fe_values_N.shape_value(j,
q_point)) *
fe_values_N.JxW(q_point));
cell_mass_matrix_N(i,j) += fe_values_N.shape_value(i
,q_point) *
fe_values_N.shape_value(j,q_point) *
fe_values_N.JxW(q_point);
cell_matrix_TE(i,j) += ((-diffusion_coefficient_TE(
fe_values_N.shape_value(i, q_point)*fe_values_N.shape_value(j, q_point)*
fe_values_N.JxW(q_point),
fe_values_TE.shape_value(i, q_point)*fe_values_TE.shape_value(j, q_point)*
fe_values_TE.JxW(q_point),
)
*
fe_values_TE.shape_grad(i,
q_point) *
fe_values_TE.shape_grad(j,
q_point)
- absorption_cross_section_TE(
fe_values_N.shape_value(i, q_point)*fe_values_N.shape_value(j, q_point)*
fe_values_N.JxW(q_point),
fe_values_TE.shape_value(i, q_point)*fe_values_TE.shape_value(j, q_point)*
fe_values_TE.JxW(q_point)) *
fe_values_TE.JxW(q_point));
cell_mass_matrix_TE(i,j) += fe_values_TE.shape_value
(i,q_point) *
fe_values_TE.shape_value(j,q_point) *
fe_values_TE.JxW(q_point);
}
cell->get_dof_indices(local_dof_indices_N);
constraint_matrix.distribute_local_to_global(cell_matrix_N,
local_dof_indices_N,system_matrix_N);
constraint_matrix.distribute_local_to_global(cell_mass_matrix_N,
local_dof_indices_N,mass_matrix_N);
cell->get_dof_indices(local_dof_indices_TL);
constraint_matrix.distribute_local_to_global(cell_matrix_TL,
local_dof_indices_TL,system_matrix_TL);
constraint_matrix.distribute_local_to_global(cell_mass_matrix_TL
,local_dof_indices_TL,mass_matrix_TL);
}
But I am not sure if that is even remotely correct. How does fe_values_N
differ from fe_values_TE? But how should I initialize them else? Should I
use FESystem instead? If yes, how should I initialize them?
Furthermore, can I use one dof_handler for both variables?
--
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.
