FEInterfaceValues<dim>

54 views
Skip to first unread message

Alberto Salvadori

unread,
Dec 16, 2020, 1:50:42 PM12/16/20
to deal.II User Group
Dear community

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.

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.

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);

Can anyone address me some reading on this notion and how it is dealt with

I appreciate your help

Alberto



Wolfgang Bangerth

unread,
Dec 16, 2020, 2:32:21 PM12/16/20
to dea...@googlegroups.com

Alberto,
have you taken a look at step-47? That's probably the program with the
"purest" use of that class.
Best
W.

On 12/16/20 11:50 AM, Alberto Salvadori wrote:
> Dear community
>
> I would like to have a few more details on the class FEInterfaceValues<dim>
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassFEInterfaceValues.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485072776%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=1Fyvgd963Qh0vJjZq5RRtEBChAZAXCa0Tjs7ayAiYZg%3D&reserved=0>.
> since it is my feeling that FEInterfaceValues<dim>
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassFEInterfaceValues.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485082802%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=wSWy88SFBciv5Qi6cwIdEbXHXTbRpe%2FwuF%2BGqUMwNjg%3D&reserved=0>
> might be conveniently used for interfaces, even out of DG methods.
>
> Specifically I'd like to figure out if this notion can be used without
> recurring to the MeshWorker::mesh_loop
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2Fgroup__MeshWorker.html%23gab75fe399736977323bb319959ab66feb&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485082802%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=AcVi0QtyMphSjCwJxLU3XjfntYNYPZt9A62%2FTdZ0OgM%3D&reserved=0>,
> which I am not quite confident at present.
>
> I'd like to play autonomously with the class FEInterfaceValues<dim>
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassFEInterfaceValues.html&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485092771%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=sO2EXLTMRqhCBJKVmtg0f%2FM%2BqxYGafBWd0LdXo4qoYE%3D&reserved=0>
> but unfortunately I cannot grasp the notion of subface, which is apparently
> required by the reinit as in step 12:
>  fe_iv.reinit
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2FclassFEInterfaceValues.html%23af6aa20b6652c2605d4cee36021f6bcd2&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485092771%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=S9u2Xv49A%2Fan0RyhBfE4EmR2%2FevhoF4coo%2BkTvt5GyU%3D&reserved=0>(cell,
> f, sf, ncell, nf, nsf);
>
> Can anyone address me some reading on this notion and how it is dealt with
> by the MeshWorker::mesh_loop
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.dealii.org%2Fcurrent%2Fdoxygen%2Fdeal.II%2Fgroup__MeshWorker.html%23gab75fe399736977323bb319959ab66feb&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485102756%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=UqH1dDvZ0NLeU9IIt1FEy1hfJAxXluP2xYkaDlQocO0%3D&reserved=0> ?
>
> I appreciate your help
>
> Alberto
>
>
>
>
>
> Informativa sulla Privacy: http://www.unibs.it/node/8155
> <https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.unibs.it%2Fnode%2F8155&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485102756%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=lWwcRecpRmaLatMfPy5PK91dgrIV%2FRjV9NjlLEZM%2F9U%3D&reserved=0>
>
> --
> The deal.II project is located at http://www.dealii.org/
> <https://nam01.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.dealii.org%2F&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485112753%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=7dOj9gukULSm9sh4zspHljWroqwpi5Kt2arkXOsaWJQ%3D&reserved=0>
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fforum%2Fdealii%3Fhl%3Den&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485112753%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=Jq8v930wtj9zxnbL68%2Fh6b9uBqcreghcSJpfGF3JnTg%3D&reserved=0>
> ---
> 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
> <mailto:dealii+un...@googlegroups.com>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/1d271fc5-1017-49a3-80bf-cdd17e2e70aan%40googlegroups.com
> <https://nam01.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fdealii%2F1d271fc5-1017-49a3-80bf-cdd17e2e70aan%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7CWolfgang.Bangerth%40colostate.edu%7C26a50c358048457e382608d8a1f37ff2%7Cafb58802ff7a4bb1ab21367ff2ecfc8b%7C0%7C0%7C637437414485122745%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=cTx3wHK7vX%2BTkiBLfYtuxLcnSILhHRlGM%2BALvgcY1kA%3D&reserved=0>.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

Timo Heister

unread,
Dec 16, 2020, 5:25:44 PM12/16/20
to dea...@googlegroups.com
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/

Alberto Salvadori

unread,
Dec 22, 2020, 4:26:22 AM12/22/20
to dea...@googlegroups.com
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



--
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.

Alberto Salvadori

unread,
Dec 23, 2020, 1:11:47 AM12/23/20
to deal.II User Group
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

Alberto Salvadori

unread,
Dec 23, 2020, 2:14:04 AM12/23/20
to deal.II User Group
Some further information:

- no errors are detected by PETSc. By using it in place of Trilinos the code seem to work OK.
- if only the rhs is distributed, no errors again. Errors are due to the matrix distribution in Trilinos. 

Thank you again
Alberto

Wolfgang Bangerth

unread,
Dec 23, 2020, 2:43:22 PM12/23/20
to dea...@googlegroups.com

Alberto,
my suspicion would be that you are computing terms for the interface matrix
that are not part of the sparsity pattern of the Trilinos matrix. (PETSc does
not require you to say up front where the nonzero entries of the matrix are
going to be, if I recall correctly.)

This can happen because for "normal" CG methods, you only get entries in the
matrix if DoFs i,j are located on the same cell. But for methods that contain
jump terms on the interfaces, you also get contributions for DoFs i,j located
on different cells (and not necessarily at the face you're integrating over).
Take a look at how step-47 builds its sparsity pattern.

Best
W.

Timo Heister

unread,
Dec 26, 2020, 5:00:22 PM12/26/20
to dea...@googlegroups.com
Alberto,

> 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

Currently, FEInterfaceValues does not have support for
get_function_values() directly but this is something I am planning to
do soon.

Your code looks reasonable, but I think you can use the existing
FEFaceValues objects inside the FEInterfaceValues like

feiv.get_fe_face_values(0).get_function_values(...)
feiv.get_fe_face_values(1).get_function_values(...)

Alberto Salvadori

unread,
Dec 27, 2020, 3:25:24 AM12/27/20
to dea...@googlegroups.com
Thank you very much, Timo.
I have my code running reasonably well now and this is mostly because of your suggestions.
The files you pointed me out to were exactly what I was looking for.
Perhaps this little effort in writing interfaces for continuous fe might worth some description in the deal.ii
library.

Thank you again
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


--
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.

Alberto Salvadori

unread,
Jan 8, 2021, 12:19:48 PM1/8/21
to dea...@googlegroups.com
Wolfgang,

thank you for the hint. I looked at setp 47 sparsity pattern. If I understood it right, 
the main difference stands in the locally_relevant_dofs, which are no longer extracted and used. 
Rather, the sparsity pattern is based on the whole dof_handler, as per the instruction

DoFTools::make_flux_sparsity_pattern( this->dof_handler, dsp, this->hanging_node_constraints, true );

At present, I am a bit confused by the final statements of the setup_system, namely:

system_matrix.reinit(sparsity_pattern);
system_rhs.reinit(dof_handler.n_dofs());

In step 47, system_matrix is of type Sparse_Matrix, whereas I am using wrappers to PETSc or TRILINOS ( i.e LA::MPI::SparseMatrix )
and they cannot be reinit with the sparsity_pattern. Similarly for system_rhs. I have been reading those classes yet could not figure out how to
circumvent the problem. Can you please address me in this regard?

Thank you so much,
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


--
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.

Wolfgang Bangerth

unread,
Jan 8, 2021, 3:10:16 PM1/8/21
to dea...@googlegroups.com
On 1/8/21 10:19 AM, Alberto Salvadori wrote:
>
> thank you for the hint. I looked at setp 47 sparsity pattern. If I understood
> it right,
> the main difference stands in the locally_relevant_dofs, which are no longer
> extracted and used.
> Rather, the sparsity pattern is based on the whole dof_handler, as per the
> instruction
>
> DoFTools::make_flux_sparsity_pattern( this->dof_handler, dsp,
> this->hanging_node_constraints, true );

No, you mistook what the difference is. The difference is that
make_sparsity_pattern() only adds entries for the matrix that correspond to
DoFs i and j that are located on the same cell.

make_flux_sparsity_pattern() also adds entries for DoFs i,j that are located
on cells that are face neighbors of each other.

In both cases, you can use the same locally_relevant_set when you set up the
sparsity pattern object.

Alberto Salvadori

unread,
Jan 9, 2021, 5:31:37 AM1/9/21
to dea...@googlegroups.com
Thank you, Wolfgang. The explanation was clear as usual, I understood the principle but overlooked the DoFTools class. This is now set.
Unfortunately, it seems that it does not solve my issue and I guess I figured out why.

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 - as depicted already (I won't repeat it here).

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 FOR TRILINOS (AND NOT FOR PETSc, which works OK).

According to your remark "make_flux_sparsity_pattern() also adds entries for DoFs i,j that are located
on cells that are face neighbors of each other" the error may arise because in defining the interfaces I am:
1 - duplicating the nodes in order to separate the initial triangulation into two
2 - building a new data structure that provide connectivity for the dofs of the separated elements, exploiting the class FEInterfaceValues
I wonder if by doing so, the notion of "cells that are face neighbors of each other" may not apply anymore and hence, as you suspected, I am computing terms for the interface matrix

that are not part of the sparsity pattern of the Trilinos matrix.

If this is correct, do you see a way to circumvent the TRILINOS problem I am facing?

Thank you,
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


--
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.

Timo Heister

unread,
Jan 9, 2021, 2:24:47 PM1/9/21
to dea...@googlegroups.com
> If this is correct, do you see a way to circumvent the TRILINOS problem I am facing?

If you need entries i,j in the sparsity pattern and none of the
make_*_sparsity_pattern() function apply, you will need to add them
manually. This should be as simple as following the same logic you are
using in the assembly, where you identify DoF indices (you called them
local_interface_dof_indices) and you add them using

dsp.add(i,j)

where i and j are all indices in local_interface_dof_indices.
Reply all
Reply to author
Forward
0 new messages