Dear Wolfgang and Daniel,
Thanks for your directions. I've reviewed Step-20, Step-60, and Step-70
and started working on my own example.
In this case, I'm trying to combine two types of dimensionality: the <3,3>
Solid Elements and the <2,3> "Shell" Elements. I've attached a picture of
the meshes. I am trying to do a block system. For the moment, I am not
trying to link the two meshes (off-diagonal blocks), only get them into the
same system matrix. Each is independently constrained with appropriate
boundary conditions. I've run the meshes separately and they work fine.
My understanding is that I cannot combine several objects because of the
difference in <dim,spacedim>. So, I've created duplicates of several
objects, using "_sld" to note for the Solid and "_sh" for the shell.
Triangulation<3> triangulation_sld;
Triangulation<2,3> triangulation_sh;
DoFHandler<3> dof_handler_sld;
DoFHandler<2,3> dof_handler_sh;
FESystem<3> fe_sld;
FESystem<2,3> fe_sh;
const MappingFE<3> mapping_sld;
const MappingFE<2,3> mapping_sh;
const QGauss<3> quadrature_sld;
const QGauss<2> quadrature_sh;
I then create the following objects that will be shared by both the Solid
and Shell elements.
AffineConstraints<double> constraints;
BlockSparsityPattern sparsity_pattern;
BlockSparseMatrix<double> system_matrix;
BlockSparseMatrix<double> preconditioner_matrix;
BlockVector<double> solution;
BlockVector<double> system_rhs;
BlockVector<double> locally_relevant_solution;
I proceed by reading two separate external meshes and attaching them to the
appropriate triangulations.
Then I begin the Setup_System(). I start by getting the number of DOFs
from each of the triangulations and adding them
in a vector of the DOFs for each block. I am saying that the Solid is
block 0 and the Shell is block 1.
const unsigned int n_dof_sld = dof_handler_sld.n_dofs(); // = 24
const unsigned int n_dof_sh = dof_handler_sh.n_dofs(); // = 48
const std::vector<unsigned int> dofs_per_block = {n_dof_sld, n_dof_sh};
I then begin work on the sparsity pattern:
BlockDynamicSparsityPattern dsp(dofs_per_block , dofs_per_block);
for (unsigned int i=0; i<dofs_per_block.size(); ++i)
for (unsigned int j=0; j<dofs_per_block.size(); ++j)
dsp.block(i,j).reinit(dofs_per_block[i], dofs_per_block[j]);
dsp.collect_sizes();
sparsity_pattern.copy_from(dsp);
system_matrix.reinit(sparsity_pattern);
system_rhs.reinit(dofs_per_block);
>From what I can see, the system_matrix looks right with n_blocks = 2 and
start_indices = {0, 24, 72}, as expected.
Next, I run through the Assemble_System() to create a local cell matrix. I
start with the Solid elements and cell_matrix_sld. I'm skipping the
calculation of the cell matrix for brevity's sake.
for (const auto &cell : dof_handler_sld.active_cell_iterators())
{
// *** Placeholder for calculating cell_matrix_sld, a 24x24 matrix ***
cell->get_dof_indices(local_dof_indices_sld);
// Add the cell matrix to the system matrix
for (unsigned int i = 0; i < dofs_per_cell_sld; ++i)
for (unsigned int j = 0; j < dofs_per_cell_sld; ++j)
system_matrix.add(local_dof_indices_sld[i],
local_dof_indices_sld[j],
cell_matrix_sld(i, j));
} // end of cell loop
The bit of code 'system_matrix.add(...)' works for when i=0 and j=0, but
fails on the next loop when j=1. I get this error:
An error occurred in line <1931> of file <.../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.
So, it appears that I'm not setting something up correctly with the
system_matrix and sparsity pattern. Step-20 uses the same
system_matrix.add(...) approach, and does not have any additional
preparation of the system_matrix that I can see.
One thing I skipped is this line from Step 20:
DoFTools::make_sparsity_pattern
<https://www.dealii.org/current/doxygen/deal.II/group__constraints.html#gaf10030e7a598def1a1fa43fb0c005a93>(dof_handler,
dsp);
I'm unable to use this line because I have the two separate DOF handlers.
Do you have any idea of what I can do here?
Best regards,
Alex
On Saturday, October 28, 2023 at 6:58:21 PM UTC-4 Wolfgang Bangerth wrote:
> On 10/27/23 12:21, Daniel Arndt wrote:
> >
> > If you can formulate your problem/approach in a way (iteratively?) that
> you
> > could solve separately on the three triangulations (using input from the
> > solution on the other triangulations), that might be easier.
>
> That corresponds to an operator splitting scheme or, in the language of
> linear
> solvers, an Uzawa iteration.
>
> The alternative is to build a block system in which each block corresponds
> to
> one triangulation (2 or 3 dimensional). You'd build each block in the same
> way
> as you would with a single triangulation, including using the same
> indexing of
> rows/columns that corresponds to the DoFHandlers on these triangulations.
> The
> tricky part is assembling the off-diagonal blocks where you'd need to work
> a
> bit harder. If you can assemble such a linear system, then the solution
> becomes both relatively easy and efficient because you can use a coupled
> solver scheme.
>
> You might also want to look at some of the "non-matching" tutorial
> programs to
> see inspirations for problems that use different triangulations at once.
>
> Best
> W.
>
> --
> ------------------------------------------------------------------------
> Wolfgang Bangerth email: bang...@colostate.edu
> www: http://www.math.colostate.edu/~bangerth/
>
>
>
--
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.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/1042bdbd-a25f-4a34-b4d5-4126c2b4523en%40googlegroups.com.
