Skip to content

Commit

Permalink
Draft of FECouplingValues.
Browse files Browse the repository at this point in the history
  • Loading branch information
@luca-heltai
luca-heltai committed Sep 6, 2023
1 parent 70cb5f8 commit 900e1ec
Show file tree
Hide file tree
Showing 13 changed files with 2,546 additions and 19 deletions.
4 changes: 4 additions & 0 deletions doc/news/changes/minor/20230721LucaHeltai
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
New: Created a new class FECouplingValues, that helps with assembling of coupled
finite element methods across different dimensions or different grids.
<br>
(Luca Heltai, 2023/07/21)
1,109 changes: 1,090 additions & 19 deletions include/deal.II/fe/fe_coupling_values.h
View file

Large diffs are not rendered by default.

4 changes: 4 additions & 0 deletions tests/fecoupling/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
cmake_minimum_required(VERSION 3.13.4)
include(../scripts/setup_testsubproject.cmake)
project(testsuite CXX)
deal_ii_pickup_tests()
224 changes: 224 additions & 0 deletions tests/fecoupling/fe_coupling_01.cc
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,224 @@
// ---------------------------------------------------------------------
//
// Copyright (C) 2023 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------


// Test basic properties of FECouplingValues

#include <deal.II/base/quadrature_lib.h>

#include <deal.II/fe/fe_coupling_values.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q.h>

#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria.h>

#include <fstream>
#include <iostream>

#include "../tests.h"

#include "../test_grids.h"


template <typename FCV>
void
inspect_fcv(const FCV & fcv,
const std::vector<types::global_dof_index> &v1,
const std::vector<types::global_dof_index> &v2)
{
deallog << "n_coupling_dofs: " << fcv.n_coupling_dofs() << std::endl;

auto indices = fcv.get_coupling_dof_indices(v1, v2);
Assert(indices.size() == fcv.n_coupling_dofs(), ExcInternalError());

deallog << "coupling_dof_indices: ";
for (const auto &i : indices)
deallog << i << ' ';
deallog << std::endl;

unsigned int idx = 0;
for (const auto &v : indices)
{
const auto pair = fcv.coupling_dof_to_dof_indices(idx);

deallog << " index " << idx << " global_dof_index: " << v
<< ", index 1: " << static_cast<int>(pair[0])
<< ", index 2: " << static_cast<int>(pair[1]) << std::endl;
++idx;
}

deallog << "n_quadrature_points: " << fcv.n_quadrature_points() << std::endl;

for (const auto &q : fcv.quadrature_point_indices())
{
const auto pair = fcv.coupling_quadrature_to_quadrature_indices(q);
const auto &[q1, q2] = fcv.quadrature_point(q);
deallog << " q index " << q << ", index 1: " << pair[0]
<< ", index 2: " << pair[1] << ", q1: " << q1 << ", q2: " << q2
<< std::endl;
}

deallog << std::endl;
}



template <int dim>
void
test()
{
Triangulation<dim> tria;
TestGrids::hyper_line(tria, 2);

DoFHandler<dim> dofh(tria);
FE_Q<dim> fe(1);
dofh.distribute_dofs(fe);

MappingQ<dim> mapping(1);
UpdateFlags update_flags = update_quadrature_points | update_JxW_values;

FEValues<dim> fv1(mapping, fe, QGauss<dim>(fe.degree + 1), update_flags);
FEValues<dim> fv2(mapping, fe, QGauss<dim>(fe.degree + 1), update_flags);

FEValues<dim> fv3(mapping, fe, QTrapezoid<dim>(), update_flags);
FEValues<dim> fv4(mapping, fe, QTrapezoid<dim>(), update_flags);

FEFaceValues<dim> fv5(mapping, fe, QTrapezoid<dim - 1>(), update_flags);
FEFaceValues<dim> fv6(mapping, fe, QTrapezoid<dim - 1>(), update_flags);

std::vector<types::global_dof_index> v1(fe.dofs_per_cell);
std::vector<types::global_dof_index> v2(fe.dofs_per_cell);

auto cell1 = dofh.begin();
fv1.reinit(cell1);
cell1->get_dof_indices(v1);

unsigned int boundary_face = numbers::invalid_unsigned_int;
unsigned int internal_face = numbers::invalid_unsigned_int;

for (const unsigned int f : cell1->face_indices())
{
if (cell1->at_boundary(f))
{
boundary_face = f;
}
else if (!cell1->at_boundary(f))
{
internal_face = f;
}
if (boundary_face != numbers::invalid_unsigned_int &&
internal_face != numbers::invalid_unsigned_int)
break;
}

auto cell2 = cell1->neighbor(internal_face);
fv2.reinit(cell2);
cell2->get_dof_indices(v2);

{
deallog << "** coupling between cell1 and cell2 (tensor) **" << std::endl;
FECouplingValues<dim> fcv(fv1, fv2, DoFCouplingType::contiguous);
inspect_fcv(fcv, v1, v2);
}
{
deallog << "** coupling between cell1 and cell2 (unrolled) **" << std::endl;
FECouplingValues<dim> fcv(fv1,
fv2,
DoFCouplingType::contiguous,
QuadratureCouplingType::unrolled);
inspect_fcv(fcv, v1, v2);
}
{
deallog << "** coupling between cell1 and cell2 Trapez (overlapping) **"
<< std::endl;
fv3.reinit(cell1);
fv4.reinit(cell2);
FECouplingValues<dim> fcv(fv3,
fv4,
DoFCouplingType::contiguous,
QuadratureCouplingType::overlapping);
inspect_fcv(fcv, v1, v2);
}
{
deallog
<< "** coupling between cell1 and boundary face, trapez (overlapping) **"
<< std::endl;
fv3.reinit(cell1);
fv5.reinit(cell1, boundary_face);
FECouplingValues<dim> fcv(fv3,
fv5,
DoFCouplingType::contiguous,
QuadratureCouplingType::overlapping);
inspect_fcv(fcv, v1, v1);
}
{
deallog
<< "** coupling between cell1 and cell1->internal_face, trapez (overlapping) **"
<< std::endl;
fv3.reinit(cell1);
fv5.reinit(cell1, internal_face);
FECouplingValues<dim> fcv(fv3,
fv5,
DoFCouplingType::contiguous,
QuadratureCouplingType::overlapping);
inspect_fcv(fcv, v1, v2);
}
{
deallog
<< "** coupling between cell1 and cell2->internal_face, trapez (overlapping) **"
<< std::endl;
fv3.reinit(cell1);
fv5.reinit(cell2, cell1->neighbor_of_neighbor(internal_face));
FECouplingValues<dim> fcv(fv3,
fv5,
DoFCouplingType::contiguous,
QuadratureCouplingType::overlapping);
inspect_fcv(fcv, v1, v2);
}
{
deallog
<< "** coupling between cell1 and cell2->other_face, trapez (overlapping) **"
<< std::endl;
fv3.reinit(cell1);
fv5.reinit(cell2, internal_face);
FECouplingValues<dim> fcv(fv3,
fv5,
DoFCouplingType::contiguous,
QuadratureCouplingType::overlapping);
inspect_fcv(fcv, v1, v2);
}
{
deallog
<< "** coupling between cell1->face(0) and cell1->face(2), trapez (tensor) **"
<< std::endl;
fv5.reinit(cell1, 0);
fv6.reinit(cell1, 2);
FECouplingValues<dim> fcv(fv5,
fv6,
DoFCouplingType::contiguous,
QuadratureCouplingType::tensor_product);
inspect_fcv(fcv, v1, v1);
}
}

int
main()
{
initlog();
test<2>();
}
137 changes: 137 additions & 0 deletions tests/fecoupling/fe_coupling_01.output
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,137 @@

DEAL::** coupling between cell1 and cell2 (tensor) **
DEAL::n_coupling_dofs: 8
DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5
DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1
DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1
DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1
DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1
DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0
DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1
DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2
DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3
DEAL::n_quadrature_points: 16
DEAL:: q index 0, index 1: 0, index 2: 0, q1: 0.211325 0.211325, q2: 1.21132 0.211325
DEAL:: q index 1, index 1: 0, index 2: 1, q1: 0.211325 0.211325, q2: 1.78868 0.211325
DEAL:: q index 2, index 1: 0, index 2: 2, q1: 0.211325 0.211325, q2: 1.21132 0.788675
DEAL:: q index 3, index 1: 0, index 2: 3, q1: 0.211325 0.211325, q2: 1.78868 0.788675
DEAL:: q index 4, index 1: 1, index 2: 0, q1: 0.788675 0.211325, q2: 1.21132 0.211325
DEAL:: q index 5, index 1: 1, index 2: 1, q1: 0.788675 0.211325, q2: 1.78868 0.211325
DEAL:: q index 6, index 1: 1, index 2: 2, q1: 0.788675 0.211325, q2: 1.21132 0.788675
DEAL:: q index 7, index 1: 1, index 2: 3, q1: 0.788675 0.211325, q2: 1.78868 0.788675
DEAL:: q index 8, index 1: 2, index 2: 0, q1: 0.211325 0.788675, q2: 1.21132 0.211325
DEAL:: q index 9, index 1: 2, index 2: 1, q1: 0.211325 0.788675, q2: 1.78868 0.211325
DEAL:: q index 10, index 1: 2, index 2: 2, q1: 0.211325 0.788675, q2: 1.21132 0.788675
DEAL:: q index 11, index 1: 2, index 2: 3, q1: 0.211325 0.788675, q2: 1.78868 0.788675
DEAL:: q index 12, index 1: 3, index 2: 0, q1: 0.788675 0.788675, q2: 1.21132 0.211325
DEAL:: q index 13, index 1: 3, index 2: 1, q1: 0.788675 0.788675, q2: 1.78868 0.211325
DEAL:: q index 14, index 1: 3, index 2: 2, q1: 0.788675 0.788675, q2: 1.21132 0.788675
DEAL:: q index 15, index 1: 3, index 2: 3, q1: 0.788675 0.788675, q2: 1.78868 0.788675
DEAL::
DEAL::** coupling between cell1 and cell2 (unrolled) **
DEAL::n_coupling_dofs: 8
DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5
DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1
DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1
DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1
DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1
DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0
DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1
DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2
DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3
DEAL::n_quadrature_points: 4
DEAL:: q index 0, index 1: 0, index 2: 0, q1: 0.211325 0.211325, q2: 1.21132 0.211325
DEAL:: q index 1, index 1: 1, index 2: 1, q1: 0.788675 0.211325, q2: 1.78868 0.211325
DEAL:: q index 2, index 1: 2, index 2: 2, q1: 0.211325 0.788675, q2: 1.21132 0.788675
DEAL:: q index 3, index 1: 3, index 2: 3, q1: 0.788675 0.788675, q2: 1.78868 0.788675
DEAL::
DEAL::** coupling between cell1 and cell2 Trapez (overlapping) **
DEAL::n_coupling_dofs: 8
DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5
DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1
DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1
DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1
DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1
DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0
DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1
DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2
DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3
DEAL::n_quadrature_points: 2
DEAL:: q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000
DEAL:: q index 1, index 1: 3, index 2: 2, q1: 1.00000 1.00000, q2: 1.00000 1.00000
DEAL::
DEAL::** coupling between cell1 and boundary face, trapez (overlapping) **
DEAL::n_coupling_dofs: 8
DEAL::coupling_dof_indices: 0 1 2 3 0 1 2 3
DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1
DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1
DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1
DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1
DEAL:: index 4 global_dof_index: 0, index 1: -1, index 2: 0
DEAL:: index 5 global_dof_index: 1, index 1: -1, index 2: 1
DEAL:: index 6 global_dof_index: 2, index 1: -1, index 2: 2
DEAL:: index 7 global_dof_index: 3, index 1: -1, index 2: 3
DEAL::n_quadrature_points: 2
DEAL:: q index 0, index 1: 0, index 2: 0, q1: 0.00000 0.00000, q2: 0.00000 0.00000
DEAL:: q index 1, index 1: 2, index 2: 1, q1: 0.00000 1.00000, q2: 0.00000 1.00000
DEAL::
DEAL::** coupling between cell1 and cell1->internal_face, trapez (overlapping) **
DEAL::n_coupling_dofs: 8
DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5
DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1
DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1
DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1
DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1
DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0
DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1
DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2
DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3
DEAL::n_quadrature_points: 2
DEAL:: q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000
DEAL:: q index 1, index 1: 3, index 2: 1, q1: 1.00000 1.00000, q2: 1.00000 1.00000
DEAL::
DEAL::** coupling between cell1 and cell2->internal_face, trapez (overlapping) **
DEAL::n_coupling_dofs: 8
DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5
DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1
DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1
DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1
DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1
DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0
DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1
DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2
DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3
DEAL::n_quadrature_points: 2
DEAL:: q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000
DEAL:: q index 1, index 1: 3, index 2: 1, q1: 1.00000 1.00000, q2: 1.00000 1.00000
DEAL::
DEAL::** coupling between cell1 and cell2->other_face, trapez (overlapping) **
DEAL::n_coupling_dofs: 8
DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5
DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1
DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1
DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1
DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1
DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0
DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1
DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2
DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3
DEAL::n_quadrature_points: 0
DEAL::
DEAL::** coupling between cell1->face(0) and cell1->face(2), trapez (tensor) **
DEAL::n_coupling_dofs: 8
DEAL::coupling_dof_indices: 0 1 2 3 0 1 2 3
DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1
DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1
DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1
DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1
DEAL:: index 4 global_dof_index: 0, index 1: -1, index 2: 0
DEAL:: index 5 global_dof_index: 1, index 1: -1, index 2: 1
DEAL:: index 6 global_dof_index: 2, index 1: -1, index 2: 2
DEAL:: index 7 global_dof_index: 3, index 1: -1, index 2: 3
DEAL::n_quadrature_points: 4
DEAL:: q index 0, index 1: 0, index 2: 0, q1: 0.00000 0.00000, q2: 0.00000 0.00000
DEAL:: q index 1, index 1: 0, index 2: 1, q1: 0.00000 0.00000, q2: 1.00000 0.00000
DEAL:: q index 2, index 1: 1, index 2: 0, q1: 0.00000 1.00000, q2: 0.00000 0.00000
DEAL:: q index 3, index 1: 1, index 2: 1, q1: 0.00000 1.00000, q2: 1.00000 0.00000
DEAL::

0 comments on commit 900e1ec

Please sign in to comment.