Dear Daniel,
Thank you very much for your rapid reply! I use overload
<https://www.dealii.org/developer/doxygen/deal.II/namespaceDoFTools.html#a9696e285b577a37522b1262e56750f3a>
as
you say, and it can run before assemble_system, the relevant part of my
code is
template <int dim>
void Problem<dim>::run ()
{
const unsigned int n_cycles = 1;
double entropy = 0;
for (unsigned int cycle=0; cycle<n_cycles; ++cycle)
{
if (cycle == 0)
{
GridGenerator::hyper_cube (triangulation_active, 0, 2.25);
triangulation_active.refine_global (7);
GridGenerator::flatten_triangulation (triangulation_active,
triangulation);
for (typename Triangulation<dim>::active_cell_iterator
cell = triangulation.begin_active();
cell != triangulation.end();
++cell)
{
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell;
++f)
{
if (cell->face(f)->at_boundary())
{
if (std::fabs(cell->face(f)->center()(0) -
(2.25)) < 1e-12)
cell->face(f)->set_boundary_id (1);
else if(std::fabs(cell->face(f)->center()(1) -
(0)) < 1e-12)
cell->face(f)->set_boundary_id (2);
else if(std::fabs(cell->face(f)->center()(1)
- (2.25)) < 1e-12)
cell->face(f)->set_boundary_id (3);
else
if(std::fabs(cell->face(f)->center()(0) - (0)) < 1e-12)
{
cell->face(f)->set_boundary_id
(0);
}
}
}
}
}
else
{
//triangulation.refine_global (1);
}
setup_dofs ();
std::cout << "start" << cycle << std::endl;
VectorTools::project (dof_handler_phi, constraints,
QGauss<dim>(degree+2),
InitialValues_phi<dim>(),
phi_0);
VectorTools::project (dof_handler_e, constraints,
QGauss<dim>(degree+2),
InitialValues_e<dim>(),
e_0);
std::cout << "phi_0" << cycle << std::endl;
//VectorTools::project (dof_handler_U, constraints_U,
QGauss<dim>(degree+2),
// InitialValues_U<dim>(),
// U_n);
//std::cout << "U_0" << cycle << std::endl;
VectorTools::project (dof_handler_q, constraints_q,
QGauss<dim>(degree+2),
InitialValues_q<dim>(),
q_0);
std::cout << "q_0" << cycle << std::endl;
assemble_system_1 (phi_0, e_0);
[...]
and my setup_dofs ():
template <int dim>
void Problem<dim>::setup_dofs ()
{
std::vector<unsigned int> block_component (4);
block_component[0] = 0;
block_component[1] = 1;
block_component[2] = 2;
block_component[3] = 3;
dof_handler.distribute_dofs (fe);
//DoFRenumbering::component_wise (dof_handler);
DoFRenumbering::component_wise (dof_handler, block_component);
constraints.clear ();
const ComponentMask component_mask = ComponentMask();
DoFTools::make_periodicity_constraints(dof_handler,0,1,0, constraints,
component_mask);
DoFTools::make_periodicity_constraints(dof_handler,2,3,1, constraints,
component_mask);
constraints.close ();
std::vector<types::global_dof_index> dofs_per_block (4);
DoFTools::count_dofs_per_block (dof_handler, dofs_per_block,
block_component);
const unsigned int n_phi = dofs_per_block[0],
n_e = dofs_per_block[1],
n_U = dofs_per_block[2],
n_q = dofs_per_block[3];
system_matrix_1.clear ();
BlockDynamicSparsityPattern dsp (4,4);
dsp.block(0,0).reinit (n_phi, n_phi);
dsp.block(1,0).reinit (n_e, n_phi);
dsp.block(2,0).reinit (n_U, n_phi);
dsp.block(3,0).reinit (n_q, n_phi);
dsp.block(0,1).reinit (n_phi, n_e);
dsp.block(1,1).reinit (n_e, n_e);
dsp.block(2,1).reinit (n_U, n_e);
dsp.block(3,1).reinit (n_q, n_e);
dsp.block(0,2).reinit (n_phi, n_U);
dsp.block(1,2).reinit (n_e, n_U);
dsp.block(2,2).reinit (n_U, n_U);
dsp.block(3,2).reinit (n_q, n_U);
dsp.block(0,3).reinit (n_phi, n_q);
dsp.block(1,3).reinit (n_e, n_q);
dsp.block(2,3).reinit (n_U, n_q);
dsp.block(3,3).reinit (n_q, n_q);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern (dof_handler, dsp, constraints, false);
sparsity_pattern.copy_from (dsp);
system_matrix_1.reinit (sparsity_pattern);
[...]
the output and error is
start0
phi_00
q_00
--------------------------------------------------------
An error occurred in line <1764> of file
</.../boost/deal.ii-8.5.1/deal.II/include/deal.II/lac/sparse_matrix.h> in
function
void
dealii::SparseMatrix<number>::add(dealii::SparseMatrix<number>::size_type,
dealii::SparseMatrix<number>::size_type, number) [with number = double;
dealii::SparseMatrix<number>::size_type = unsigned int]
The violated condition was:
(index != SparsityPattern::invalid_entry) || (value == number())
Additional information:
You are trying to access the matrix entry with index <0,1>, but this
entry does not exist in the sparsity pattern of this matrix.
The most common cause for this problem is that you used a method to build
the sparsity pattern that did not (completely) take into account all of the
entries you will later try to write into. An example would be building a
sparsity pattern that does not include the entries you will write into due
to constraints on degrees of freedom such as hanging nodes or periodic
boundary conditions. In such cases, building the sparsity pattern will
succeed, but you will get errors such as the current one at one point or
other when trying to write into the entries of the matrix.
I wander how to build the sparsity pattern completely (when I use the same
part as step-45, i.e. the same setup-dof(), but receive the same error as
above), maybe I need to add some information to the "
BlockDynamicSparsityPattern
dsp", but how to do it? I can't find a correct example to describe this.
Thank you very much!
Best,
Chucui
--
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.
