Dear community
After Timo's suggestions I attempted at assembling matrix and rhs for
problems with an interface, using continuum elements (not DG).
I mean at this: take two domains connected by springs and pull them apart,
so to separate the domains and elongate the spring.
Since the forces provided by the springs are due to the separation of the
faces, there is a contribution in the Newton Raphson system
provided by the interface.
The code I wrote is actually "almost" working, some help in addressing the
hopefully last piece of the puzzle is very much appreciated. In brief, I am
building an
interface matrix and an interface rhs in the very same way cell matrix and
cell rhs have been built in step 18:
FEInterfaceValues<dim> fiv(fe, face_quadrature_formula,
update_values | update_quadrature_points |
update_normal_vectors | update_JxW_values);
reinit_FeInterfaceValues ( fiv );
unsigned dofs_per_interfaces = fiv.n_current_interface_dofs() ;
std::vector<types::global_dof_index> local_interface_dof_indices
(dofs_per_interfaces);
FullMatrix<double> interface_matrix ( dofs_per_interfaces,
dofs_per_interfaces );
Vector<double> interface_rhs ( dofs_per_interfaces );
interface_matrix = 0.0;
interface_rhs = 0.0;
for (unsigned int q_point=0; q_point < q_points.size(); ++q_point)
{
for (unsigned int ii=0; ii<dofs_per_interfaces; ++ii)
{
for ( unsigned int jj=0; jj<dofs_per_interfaces; ++jj )
interface_matrix( ii,jj ) += ...
interface_rhs(ii) += ...
}
}
so far, this seems to work OK (although so far I did not checked if
processors own cells). Issues arise when I attempt at distributing
the matrices into system matrix: I coded this
local_interface_dof_indices = fiv.get_interface_dof_indices() ;
this->hanging_node_constraints.distribute_local_to_global(interface_matrix,
interface_rhs, local_interface_dof_indices, this->system_matrix,
this->system_rhs);
using the very same hanging_node_constraints I used for the usual continuum
elements. This last statement causes an error, though. Either running in
debug or release mode, these are the info from deal.ii:
* ----------------------------------------------------*
*Exception on processing: *
*--------------------------------------------------------*
*An error occurred in line <1683> of file
</var/folders/8z/hlb6vc015qjggytkxn84m6_c0000gn/T/heltai/spack-stage/spack-stage-dealii-9.2.0-yv7ehpigjxtbrhqi7y4menbu532z6rxi/spack-src/source/lac/trilinos_sparse_matrix.cc>
in function*
* void dealii::TrilinosWrappers::SparseMatrix::add(const
dealii::TrilinosWrappers::SparseMatrix::size_type, const
dealii::TrilinosWrappers::SparseMatrix::size_type, const
dealii::TrilinosWrappers::SparseMatrix::size_type *, const
dealii::TrilinosScalar *, const bool, const bool)*
*The violated condition was: *
* ierr == 0*
*Additional information: *
* An error with error number 2 occurred while calling a Trilinos
function*
*--------------------------------------------------------*
*Aborting!*
*----------------------------------------------------*
The hanging_node_constraints was defined in the setup_system as:
this->hanging_node_constraints.clear ();
DoFTools::make_hanging_node_constraints (this->dof_handler,
this->hanging_node_constraints);
this->hanging_node_constraints.close ();
Can I ask from some hints?
Thank you so much
Alberto
Il giorno martedì 22 dicembre 2020 alle 10:26:22 UTC+1 Alberto Salvadori ha
scritto:
> Timo,
> thank you. These two files are extremely helpful and I guess I understood
> a lot from them. Extremely well written, by the way.
> I really appreciate it.
>
> I would like to associate to each gauss point on the interface,
> identified by the command
> const auto &q_points = fiv.get_quadrature_points();
> the values of the unknown fields at the two faces of the two cells that
> form the interface itself.
> It should not be too difficult I think.
> Perhaps I should use
>
> FEFaceValues<dim> fe_face0_values (fe, face_quadrature_formula,
> update_values |
> update_quadrature_points |
> update_normal_vectors |
> update_JxW_values);
> FEFaceValues<dim> fe_face1_values (fe, face_quadrature_formula,
> update_values |
> update_quadrature_points |
> update_normal_vectors |
> update_JxW_values);
>
> fe_face0_values.reinit (cell0, face_number_on_cell0);
> fe_face1_values.reinit (cell0, face_number_on_cell1);
>
> std::vector< Vector< double > > old_solution_values_at_face_0 (
> face_quadrature_formula.size(), Vector< double >(dofs.GPDofs) );
> fe_face0_values.get_function_values ( accumulated_displacement,
> old_solution_values_at_face_0 );
> std::vector< Vector< double > > old_solution_values_at_face_1 (
> face_quadrature_formula.size(), Vector< double >(dofs.GPDofs) );
> fe_face1_values.get_function_values ( accumulated_displacement,
> old_solution_values_at_face_0 );
>
> I am concerned about the correspondence of quadrature points on the
> interface and on the cell faces. Is there a simple and effective way to
> find their association?
>
> Thanks,
> Alberto
>
>
> *Alberto Salvadori* Dipartimento di Ingegneria Meccanica e Industriale
> (DIMI)
> Universita` di Brescia, via Branze 43, 25123 Brescia
> Italy
> tel 030 3715426
>
> e-mail:
> alberto....@unibs.it
> web-page:
> http://m4lab.unibs.it/faculty.html
>
>
>
> On Wed, Dec 16, 2020 at 11:25 PM Timo Heister <hei...@clemson.edu> wrote:
>
>> Alberto,
>>
>> > I would like to have a few more details on the class
>> FEInterfaceValues<dim>.
>> > since it is my feeling that FEInterfaceValues<dim>
>> > might be conveniently used for interfaces, even out of DG methods.
>>
>> Yes, I would think so.
>>
>> > Specifically I'd like to figure out if this notion can be used without
>> recurring to the MeshWorker::mesh_loop,
>> > which I am not quite confident at present.
>>
>> Note that we are not using the MeshWorker machinery but a single
>> function that contains the loops over the cells. The code of the
>> function is "relatively" easy to understand and you will appreciate
>> the logic for the following things you would have to deal with when
>> looping manually:
>> - assembling from the finer cell for interfaces with adaptive refinement
>> - correctly handling each facet once (or twice if desired) even in
>> parallel
>> - handling of periodic boundaries
>>
>> > I'd like to play autonomously with the class FEInterfaceValues<dim>
>> > but unfortunately I cannot grasp the notion of subface, which is
>> apparently
>> > required by the reinit as in step 12:
>> > fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
>>
>> With adaptivity you will need to pass the subface as returned by
>> neighbor_of_coarser_neighbor, otherwise it is
>> numbers::invalid_unsigned_int. See
>>
>> https://github.com/dealii/dealii/blob/691ff4907964605dd21e94f06c880786af2cd612/tests/feinterface/fe_interface_values_03.cc#L94-L100
>> and
>>
>> https://github.com/dealii/dealii/blob/691ff4907964605dd21e94f06c880786af2cd612/tests/feinterface/fe_interface_values_02.cc#L88-L93
>>
>> > Can anyone address me some reading on this notion and how it is dealt
>> with
>> > by the MeshWorker::mesh_loop ?
>>
>> Take a look at the tests in tests/feinterface/:
>> https://github.com/dealii/dealii/tree/master/tests/feinterface
>>
>> Many of them use FEInterfaceValues directly on a small test mesh.
>>
>> --
>> Timo Heister
>> http://www.math.clemson.edu/~heister/
>>
>> --
>> 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+un...@googlegroups.com.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/CAMRj59Fq2wUezf%3DikmWF5_BPB-Q1zzYgjOgd4XESUYdwX_pJ_A%40mail.gmail.com
>> .
>>
>
>
> Informativa sulla Privacy: http://www.unibs.it/node/8155
>
--
Informativa sulla Privacy: http://www.unibs.it/node/8155
<http://www.unibs.it/node/8155>
--
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/95882a87-b1e3-474f-9b09-26b37a68843en%40googlegroups.com.
