Dear Martin,
after one week on a conference and two additional weeks on holidays I am
finally back in office and realised the test case. I attached the file to
the post and also log files for runs for both in Debug and Release mode.
You can call DynamicSparsityPattern::n_nonzero_elements() to get the number
> of nonzero entries in the dynamic sparsity pattern. This method also exists
> in BlockSparsityPattern (or all sparsity patterns that inherit from
> BlockSparsityPatternBase):
>
> https://dealii.org/developer/doxygen/deal.II/classBlockSparsityPatternBase.html
>
> What I'm trying to understand here is what kind of properties your problem
> has - whether there are many nonzero entries per row and other special
> things that could explain your problems.
>
> I just checked the 3D case of step-22 for the performance of
> BlockSparsityPattern::copy_from(BlockDynamicSparsityPattern) and the
> performance looks where I would expect it to be. It takes 1.19s to copy the
> sparsity pattern for a case with 1.6m DoFs (I have some modifications for
> the mesh compared to what you find online) on my laptop. Given that there
> are 275m nonzero entries in that matrix and I need to touch around 4.4 GB
> (= 4 x 275m x 4 bytes per unsigned int, once for clearing the data in the
> pattern, once for reading in the dynamic pattern, once for writing into the
> fixed pattern plus once for write-allocate on that last operation) of
> memory here, I reach 26% of the theoretical possible on this machine (~14
> GB/s memory transfer per core). While I would know how to reach more than
> 80% of peak memory bandwidth here, this function is no way near being
> relevant in the global run time in any of my performance profiles. And I'm
> likely the deal.II person with most affinity to performance numbers.
>
>
I checked the number of nonzero elements. For 133,098 dofs I got 6,347,893
nonzero elements in the BlockDynamicSparsityPattern and it took me 20.1535s
in Release and 56.5698s in Debug mode respectively to copy them.
> Is there anything else special about your configuration or problem as
> compared to the cases presented in the tutorial? What deal.II version are
> you using, what is the finite element? Any special constraints on those
> systems?
>
The configuration is special in my case. I have got two triangulations to
realize a domain which has more than one hanging node on a face. So I have
two triangulations describing the same geometry and the geometry is
separated in two parts. On the one part I have a quiet coarse mesh and on
the other part I have fine mesh. With the first triangulation I use
trilinear finite elements on the coarse part an FENothing on the finer part
and with the other triangulation I do the same thing vice versa. On the
interface I use Mortar basis functions to connect both parts in a weak
sense. Therefore I have to figure out the interface part first, then find
the corresponding degrees of freedom on the interface and with these
information I create the sparsity patterns for the connection blocks on my
own. This is why my block matrix has 3x3 blocks, the ones belonging to the
first and second triangulation and the connecting one belonging to the
Mortar part.
This is a rough overview concerning the setup situation. If we come to the
conclusion, that this is an important fact, I can explain the whole idea in
more detail. But this will be a lot of lines to write then.
Best regards
Dustin
--
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.
For more options, visit https://groups.google.com/d/optout.
#include <deal.II/base/logstream.h>
#include <deal.II/base/timer.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/hp/fe_values.h>
#include <deal.II/hp/q_collection.h>
#include <deal.II/lac/block_sparsity_pattern.h>
#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/numerics/vector_tools.h>
#include <fstream>
struct FacePointData
{
unsigned int
master_face_no,
slave_face_no,
master_q_point_index;
dealii::hp::DoFHandler<3>::cell_iterator
master_cell,
slave_cell;
std::vector<dealii::types::global_dof_index>
slave_face_dofs,
master_face_dofs,
slave_local_face_dofs,
master_local_face_dofs,
biorthogonal_face_dofs;
};
class InterfaceData
{
public:
InterfaceData ();
std::map<dealii::types::global_dof_index, dealii::IndexSet>&
biorthogonal_master_dof_connections ();
std::map<dealii::types::global_dof_index, dealii::IndexSet>&
biorthogonal_slave_dof_connections ();
void
initialise (unsigned int const n_interface_q_points);
std::map<dealii::types::global_dof_index, dealii::IndexSet>&
master_biorthogonal_dof_connections ();
std::vector<dealii::Point<3>>&
master_q_points ();
std::map<dealii::types::global_dof_index, dealii::IndexSet>&
master_slave_dof_connections ();
unsigned int
n_cells_master ();
unsigned int
n_cells_slave ();
unsigned int
ndofs_biorthogonal ();
std::vector<FacePointData>&
q_point_data ();
void
setup_data (dealii::hp::DoFHandler<3> const& dofhandler_m,
dealii::hp::DoFHandler<3> const& dofhandler_s,
dealii::hp::QCollection<2> const& face_q_collection);
std::map<dealii::types::global_dof_index, dealii::IndexSet>&
slave_biorthogonal_dof_connections ();
std::map<dealii::types::global_dof_index, dealii::IndexSet>&
slave_master_dof_connections ();
private:
void
add_interface_q_point_data (unsigned int const& index,
unsigned int const& slave_face_no,
dealii::hp::DoFHandler<3>::cell_iterator const& master_cell,
dealii::hp::DoFHandler<3>::cell_iterator const& slave_cell);
void
detect_contact_interface(dealii::hp::DoFHandler<3> const& dofhandler_m,
dealii::hp::DoFHandler<3> const& dofhandler_s,
dealii::hp::QCollection<2> const& face_q_collection,
std::vector<dealii::hp::DoFHandler<3>::active_cell_iterator>& master_cells,
std::vector<dealii::hp::DoFHandler<3>::active_cell_iterator>& slave_cells,
std::vector<unsigned int>& slave_face_no,
unsigned int& master_counter,
unsigned int& slave_counter,
unsigned int& n_interface_q_points);
void
setup_interface_q_point_data(std::vector<dealii::hp::DoFHandler<3>::active_cell_iterator> const& master_cells,
std::vector<dealii::hp::DoFHandler<3>::active_cell_iterator> const& slave_cells,
std::vector<unsigned int> const& slave_face_no,
unsigned int const n_master_cells,
unsigned int const n_slave_cells,
dealii::hp::FECollection<3> const& fe_collection,
dealii::hp::QCollection<2> const& face_q_collection);
void
setup_interface_dof_connections (unsigned int const ndofs_m,
unsigned int const ndofs_s);
unsigned int
ndofs_bi,
ncells_m,
ncells_s;
std::vector<dealii::Point<3>>
m_q_points;
std::vector<FacePointData>
qpoint_data;
std::map<dealii::types::global_dof_index, dealii::IndexSet>
m_bi_dof_conn,
s_bi_dof_conn,
bi_m_dof_conn,
bi_s_dof_conn,
m_s_dof_conn,
s_m_dof_conn;
};
using namespace dealii;
InterfaceData::InterfaceData ()
:
ndofs_bi(0),
ncells_m (0),
ncells_s (0)
{
}
std::map<types::global_dof_index, IndexSet>&
InterfaceData::biorthogonal_master_dof_connections ()
{
return bi_m_dof_conn;
}
std::map<types::global_dof_index, IndexSet>&
InterfaceData::biorthogonal_slave_dof_connections ()
{
return bi_s_dof_conn;
}
void
InterfaceData::initialise (unsigned int const n_interface_q_points)
{
{
std::vector<FacePointData> tmp;
tmp.swap(qpoint_data);
qpoint_data.resize(n_interface_q_points);
}
{
std::vector<Point<3>> tmp;
tmp.swap(m_q_points);
m_q_points.resize(n_interface_q_points);
}
std::map<types::global_dof_index, IndexSet>
tmp_m_bi,
tmp_s_bi,
tmp_bi_m,
tmp_bi_s,
tmp_m_s,
tmp_s_m;
tmp_m_bi.swap(m_bi_dof_conn);
tmp_s_bi.swap(s_bi_dof_conn);
tmp_bi_m.swap(bi_m_dof_conn);
tmp_bi_s.swap(bi_s_dof_conn);
tmp_m_s.swap(m_s_dof_conn);
tmp_s_m.swap(s_m_dof_conn);
}
std::map<types::global_dof_index, IndexSet>&
InterfaceData::master_biorthogonal_dof_connections ()
{
return m_bi_dof_conn;
}
std::vector<Point<3>>&
InterfaceData::master_q_points()
{
Assert (m_q_points.size() > 0,
ExcNotInitialized ());
return m_q_points;
}
std::map<types::global_dof_index, IndexSet>&
InterfaceData::master_slave_dof_connections ()
{
return m_s_dof_conn;
}
unsigned int
InterfaceData::n_cells_master ()
{
return ncells_m;
}
unsigned int
InterfaceData::n_cells_slave ()
{
return ncells_s;
}
unsigned int
InterfaceData::ndofs_biorthogonal()
{
Assert (ndofs_bi > 0,
ExcInternalError ("No degrees of freedom found on the interface."));
return ndofs_bi;
}
std::vector<FacePointData>&
InterfaceData::q_point_data()
{
Assert(qpoint_data.size() > 0,
ExcNotInitialized());
return qpoint_data;
}
void
InterfaceData::setup_data (hp::DoFHandler<3> const& dofhandler_m,
hp::DoFHandler<3> const& dofhandler_s,
hp::QCollection<2> const& face_q_collection)
{
std::vector<hp::DoFHandler<3>::active_cell_iterator>
master_cells(dofhandler_m.get_triangulation().n_active_cells(),
dofhandler_m.end()),
slave_cells (dofhandler_s.get_triangulation().n_active_cells (),
dofhandler_s.end ());
std::vector<unsigned int>
slave_face_no (dofhandler_s.get_triangulation().n_faces(),
7);
unsigned int
slave_counter = 0,
master_counter = 0,
n_interface_q_points = 0;
detect_contact_interface(dofhandler_m,
dofhandler_s,
face_q_collection,
master_cells,
slave_cells,
slave_face_no,
master_counter,
slave_counter,
n_interface_q_points);
initialise(n_interface_q_points);
hp::FECollection<3> const&
fe_collection = dofhandler_m.get_fe();
setup_interface_q_point_data(master_cells,
slave_cells,
slave_face_no,
master_counter,
slave_counter,
fe_collection,
face_q_collection);
setup_interface_dof_connections(dofhandler_m.n_dofs(),
dofhandler_s.n_dofs());
}
std::map<types::global_dof_index, IndexSet>&
InterfaceData::slave_biorthogonal_dof_connections ()
{
return s_bi_dof_conn;
}
std::map<types::global_dof_index, IndexSet>&
InterfaceData::slave_master_dof_connections ()
{
return s_m_dof_conn;
}
void
InterfaceData::add_interface_q_point_data (unsigned int const& index,
unsigned int const& slave_face_no,
dealii::hp::DoFHandler<3>::cell_iterator const& master_cell,
dealii::hp::DoFHandler<3>::cell_iterator const& slave_cell)
{
FiniteElement<3> const*
master_fe(&master_cell->get_fe());
FiniteElement<3> const*
slave_fe(&slave_cell->get_fe());
unsigned int const
master_face_no = 1 - 2*(slave_face_no % 2) + slave_face_no;
qpoint_data[index].master_q_point_index = index;
qpoint_data[index].master_face_no = master_face_no;
qpoint_data[index].master_cell = master_cell;
qpoint_data[index].slave_face_no = slave_face_no;
qpoint_data[index].slave_cell = slave_cell;
qpoint_data[index].master_face_dofs.resize(master_fe->dofs_per_face);
qpoint_data[index].slave_face_dofs.resize(slave_fe->dofs_per_face);
master_cell->face(master_face_no)->get_dof_indices(qpoint_data[index].master_face_dofs,
master_cell->active_fe_index());
slave_cell->face(slave_face_no)->get_dof_indices(qpoint_data[index].slave_face_dofs,
slave_cell->active_fe_index());
for (unsigned int i = 0;
i < master_fe->dofs_per_cell;
++i)
{
if (master_fe->has_support_on_face(i,master_face_no))
{
qpoint_data[index].master_local_face_dofs.push_back(i);
}
}
for (unsigned int i = 0;
i < slave_fe->dofs_per_cell;
++i)
{
if (slave_fe->has_support_on_face(i, slave_face_no))
{
qpoint_data[index].slave_local_face_dofs.push_back(i);
}
}
}
void
InterfaceData::detect_contact_interface(dealii::hp::DoFHandler<3> const& dofhandler_m,
dealii::hp::DoFHandler<3> const& dofhandler_s,
dealii::hp::QCollection<2> const& face_q_collection,
std::vector<dealii::hp::DoFHandler<3>::active_cell_iterator>& master_cells,
std::vector<dealii::hp::DoFHandler<3>::active_cell_iterator>& slave_cells,
std::vector<unsigned int>& slave_face_no,
unsigned int& master_counter,
unsigned int& slave_counter,
unsigned int& n_interface_q_points)
{
ncells_m = 0;
ncells_s = 0;
typename hp::DoFHandler<3>::active_cell_iterator
cell = dofhandler_s.begin_active(),
endc = dofhandler_s.end();
for (;
cell != endc;
++cell)
{
if (dynamic_cast<const FE_Nothing<3>*>(&cell->get_fe().base_element(0)))
{
continue;
}
cell->set_user_flag();
++ncells_s;
for (unsigned int face_no = 0;
face_no < GeometryInfo<3>::faces_per_cell;
++face_no)
{
if (cell->face(face_no)->at_boundary())
{
continue;
}
if (dynamic_cast<const FE_Nothing<3>*>(&cell->neighbor(face_no)->get_fe().base_element(0)))
{
slave_cells[slave_counter] = cell;
slave_face_no[slave_counter] = face_no;
n_interface_q_points += face_q_collection[cell->active_fe_index()].size();
++slave_counter;
}
}
}
cell = dofhandler_m.begin_active(),
endc = dofhandler_m.end();
for (;
cell != endc;
++cell)
{
if (dynamic_cast<const FE_Nothing<3>*>(&cell->get_fe().base_element(0)))
{
continue;
}
cell->set_user_flag();
++ncells_m;
for (unsigned int face_no = 0;
face_no < GeometryInfo<3>::faces_per_cell;
++face_no)
{
if (cell->face(face_no)->at_boundary())
{
continue;
}
if (dynamic_cast<const FE_Nothing<3>*>(&cell->neighbor(face_no)->get_fe().base_element(0)))
{
master_cells[master_counter] = cell;
++master_counter;
}
}
}
}
void
InterfaceData::setup_interface_q_point_data(std::vector<dealii::hp::DoFHandler<3>::active_cell_iterator> const& master_cells,
std::vector<dealii::hp::DoFHandler<3>::active_cell_iterator> const& slave_cells,
std::vector<unsigned int> const& slave_face_no,
unsigned int const n_master_cells,
unsigned int const n_slave_cells,
dealii::hp::FECollection<3> const& fe_collection,
dealii::hp::QCollection<2> const& face_q_collection)
{
unsigned int
index = 0;
hp::FEFaceValues<3> fe_face_values (fe_collection,
face_q_collection,
update_quadrature_points);
for (unsigned int i = 0;
i < n_slave_cells;
++i)
{
slave_cells[i]->face(slave_face_no[i])->set_user_pointer(&qpoint_data[index]);
fe_face_values.reinit (slave_cells[i],
slave_face_no[i]);
std::vector<Point<3>>
q_points(fe_face_values.get_present_fe_values().get_quadrature_points());
for (unsigned int q = 0;
q < q_points.size();
++q)
{
for (unsigned int j = 0;
j < n_master_cells;
++j)
{
if (master_cells[j]->point_inside(q_points[q]))
{
add_interface_q_point_data(index,
slave_face_no[i],
master_cells[j],
slave_cells[i]);
m_q_points[index] = StaticMappingQ1<3>::mapping.transform_real_to_unit_cell (master_cells[j],
q_points[q]);
}
}
++index;
}
}
Assert(index == qpoint_data.size (),
ExcInternalError());
}
void
InterfaceData::setup_interface_dof_connections (unsigned int const ndofs_m,
unsigned int const ndofs_s)
{
std::map<types::global_dof_index, std::pair<types::global_dof_index, IndexSet>>
s_s_dof_conn;
unsigned int
counter = 0;
for (unsigned int i = 0;
i < qpoint_data.size();
++i)
{
std::vector<types::global_dof_index>
master_dofs (qpoint_data[i].master_face_dofs),
slave_dofs (qpoint_data[i].slave_face_dofs);
for (unsigned int k = 0;
k < master_dofs.size();
++k)
{
for (unsigned int l = 0;
l < slave_dofs.size();
++l)
{
if (m_s_dof_conn.find(master_dofs[k]) == m_s_dof_conn.end())
{
m_s_dof_conn.insert(std::make_pair(master_dofs[k],
IndexSet (ndofs_s)));
m_s_dof_conn[master_dofs[k]].add_index(slave_dofs[l]);
}
else
{
m_s_dof_conn[master_dofs[k]].add_index(slave_dofs[l]);
m_s_dof_conn[master_dofs[k]].compress();
}
if (s_m_dof_conn.find(slave_dofs[l]) == s_m_dof_conn.end())
{
s_m_dof_conn.insert(std::make_pair(slave_dofs[l],
IndexSet (ndofs_m)));
s_m_dof_conn[slave_dofs[l]].add_index(master_dofs[k]);
s_s_dof_conn.insert(std::make_pair(slave_dofs[l],
std::make_pair(counter, IndexSet (ndofs_s))));
s_s_dof_conn[slave_dofs[l]].second.add_indices(slave_dofs.begin(),
slave_dofs.end());
++counter;
}
else
{
s_m_dof_conn[slave_dofs[l]].add_index(master_dofs[k]);
s_m_dof_conn[slave_dofs[k]].compress();
s_s_dof_conn[slave_dofs[l]].second.add_indices(slave_dofs.begin(),
slave_dofs.end());
s_s_dof_conn[slave_dofs[k]].second.compress();
}
}
}
}
ndofs_bi = s_s_dof_conn.size();
std::map<types::global_dof_index, std::pair<types::global_dof_index, IndexSet>>::iterator
iter = s_s_dof_conn.begin();
for (;
iter != s_s_dof_conn.end();
++iter)
{
bi_s_dof_conn[iter->second.first] = iter->second.second;
bi_s_dof_conn[iter->second.first].compress();
bi_m_dof_conn[iter->second.first] = s_m_dof_conn[iter->first];
bi_m_dof_conn[iter->second.first].compress();
}
std::map<types::global_dof_index, IndexSet>::iterator
it;
for (it = bi_s_dof_conn.begin();
it != bi_s_dof_conn.end();
++it)
{
for (unsigned int i = 0;
i < it->second.n_elements();
++i)
{
unsigned int const
s_dof = it->second.nth_index_in_set(i);
if (s_bi_dof_conn.find(s_dof) == s_bi_dof_conn.end())
{
s_bi_dof_conn.insert(std::make_pair(s_dof,
IndexSet(ndofs_bi)));
s_bi_dof_conn[s_dof].add_index(it->first);
}
else
{
s_bi_dof_conn[s_dof].add_index(it->first);
s_bi_dof_conn[s_dof].compress();
}
}
}
for (it = bi_m_dof_conn.begin();
it != bi_m_dof_conn.end();
++it)
{
for (unsigned int i = 0;
i < it->second.n_elements();
++i)
{
unsigned int const
m_dof = it->second.nth_index_in_set(i);
if (m_bi_dof_conn.find(m_dof) == m_bi_dof_conn.end())
{
m_bi_dof_conn.insert(std::make_pair(m_dof,
IndexSet(ndofs_bi)));
m_bi_dof_conn[m_dof].add_index(it->first);
}
else
{
m_bi_dof_conn[m_dof].add_index(it->first);
m_bi_dof_conn[m_dof].compress();
}
}
}
for (unsigned int i = 0;
i < qpoint_data.size();
++i)
{
unsigned int const
n_face_dofs = qpoint_data[i].slave_face_dofs.size();
qpoint_data[i].biorthogonal_face_dofs.resize(n_face_dofs);
Assert (n_face_dofs > 0,
ExcNotInitialized ());
for (unsigned int j = 0;
j < n_face_dofs;
++j)
{
unsigned int const
s_dof = qpoint_data[i].slave_face_dofs[j];
qpoint_data[i].biorthogonal_face_dofs[j] = s_s_dof_conn[s_dof].first;
}
}
}
class HybridMesh
{
public:
HybridMesh ();
~HybridMesh ();
void
run ();
private:
void
make_biorthogonal_sparsity_pattern (dealii::BlockDynamicSparsityPattern& dsp);
void
make_grid ();
void
setup_intial_fe_distribution();
void
setup_system ();
dealii::Triangulation<3>
triangulation_m,
triangulation_s;
dealii::hp::DoFHandler<3>
dofhandler_m,
dofhandler_s;
dealii::hp::FECollection<3>
fe_collection;
dealii::hp::QCollection<3>
q_collection;
dealii::hp::QCollection<2>
face_q_collection;
dealii::BlockSparsityPattern
sparsity_pattern;
dealii::BlockSparseMatrix<double>
system_matrix;
dealii::ConstraintMatrix
constraints,
constraints_m,
constraints_s;
InterfaceData interface_data;
std::vector<FacePointData>
data_interface;
dealii::TimerOutput
timer;
unsigned int const
tdim;
unsigned int
ndofs_m,
ndofs_s,
n_cells_m,
n_cells_s;
};
HybridMesh::HybridMesh ()
:
dofhandler_m (triangulation_m),
dofhandler_s (triangulation_s),
q_collection (QGauss<3> (2)),
face_q_collection (QGauss<2> (2)),
timer (deallog.get_file_stream(),
TimerOutput::every_call_and_summary,
TimerOutput::wall_times),
tdim (3),
ndofs_m (0),
ndofs_s (0),
n_cells_m (0),
n_cells_s (0)
{
fe_collection.push_back(FESystem<3> (FE_Q<3> (1), tdim));
fe_collection.push_back(FESystem<3> (FE_Nothing<3> (1), tdim));
}
HybridMesh::~HybridMesh()
{
dofhandler_m.clear();
dofhandler_s.clear();
}
void
HybridMesh::run()
{
make_grid();
timer.enter_section("Setup");
timer.enter_subsection("setup_intial_fe_distribution");
setup_intial_fe_distribution();
timer.leave_subsection("setup_intial_fe_distribution");
timer.enter_subsection("setup_system");
setup_system();
timer.leave_subsection("setup_system");
}
void
HybridMesh::make_biorthogonal_sparsity_pattern (BlockDynamicSparsityPattern& dsp)
{
std::map<types::global_dof_index, IndexSet>::iterator
it;
for (it = interface_data.master_biorthogonal_dof_connections().begin();
it != interface_data.master_biorthogonal_dof_connections().end();
++it)
{
std::vector<unsigned int>
entries;
it->second.fill_index_vector(entries);
dsp.block(0,2).add_entries(it->first,
entries.begin(),
entries.end(),
true);
}
for (it = interface_data.slave_biorthogonal_dof_connections().begin();
it != interface_data.slave_biorthogonal_dof_connections().end();
++it)
{
std::vector<unsigned int>
entries;
it->second.fill_index_vector(entries);
dsp.block(1,2).add_entries(it->first,
entries.begin(),
entries.end(),
true);
}
for (it = interface_data.biorthogonal_master_dof_connections().begin();
it != interface_data.biorthogonal_master_dof_connections().end();
++it)
{
std::vector<unsigned int>
entries;
it->second.fill_index_vector(entries);
dsp.block(2,0).add_entries(it->first,
entries.begin(),
entries.end(),
true);
}
for (it = interface_data.biorthogonal_slave_dof_connections().begin();
it != interface_data.biorthogonal_slave_dof_connections().end();
++it)
{
std::vector<unsigned int>
entries;
it->second.fill_index_vector(entries);
dsp.block(2,1).add_entries(it->first,
entries.begin(),
entries.end(),
true);
}
}
void
HybridMesh::make_grid()
{
Point<3> const
p1 (1.0, 1.0, 1.0),
p2 (11.0, 4.0, 1.05),
p3 (1.0, 1.0, 1.0),
p4 (11.0, 4.0, 1.05);
std::vector<unsigned int>
rep1 = {10,4,1},
rep2 = {200,60,1};
GridGenerator::subdivided_hyper_rectangle(triangulation_m,
rep1,
p1,
p2);
GridGenerator::subdivided_hyper_rectangle(triangulation_s,
rep2,
p1,
p2);
typename Triangulation<3>::active_cell_iterator
cell = triangulation_m.begin_active(),
endc = triangulation_m.end();
for (;
cell != endc;
++cell)
{
for (unsigned int face_no = 0;
face_no < GeometryInfo<3>::faces_per_cell;
++face_no)
{
if (cell->face(face_no)->center()[0] == 11.0)
cell->face(face_no)->set_boundary_id(2);
if ((cell->face(face_no)->center()[2] == 1.05) &&
(cell->face(face_no)->center()[0] > 6.0))
cell->face(face_no)->set_boundary_id(4);
}
}
cell = triangulation_s.begin_active();
endc = triangulation_s.end();
for (;
cell != endc;
++cell)
{
for (unsigned int face_no = 0;
face_no < GeometryInfo<3>::faces_per_cell;
++face_no)
{
if (cell->face(face_no)->center()[0] == 1.0)
cell->face(face_no)->set_boundary_id(1);
if ((cell->face(face_no)->center()[2] == 1.0) &&
(cell->face(face_no)->center()[0] < 3.0))
cell->face(face_no)->set_boundary_id(3);
if ((cell->face(face_no)->center()[2] == 1.05) &&
(cell->face(face_no)->center()[0] < 3.0))
cell->face(face_no)->set_boundary_id(3);
}
}
triangulation_s.refine_global();
}
void
HybridMesh::setup_intial_fe_distribution()
{
typename hp::DoFHandler<3>::active_cell_iterator
cell = dofhandler_m.begin_active(),
endc = dofhandler_m.end ();
for (;
cell != endc;
++cell)
{
if (cell->center()[0] < 4.0)
{
cell->set_active_fe_index(1);
}
else
{
cell->set_active_fe_index(0);
}
}
cell = dofhandler_s.begin_active(),
endc = dofhandler_s.end ();
for (;
cell != endc;
++cell)
{
if (cell->center()[0] < 4.0)
{
cell->set_active_fe_index(0);
}
else
{
cell->set_active_fe_index(1);
}
}
}
void
HybridMesh::setup_system()
{
system_matrix.clear();
dofhandler_m.distribute_dofs(fe_collection);
dofhandler_s.distribute_dofs(fe_collection);
interface_data.setup_data (dofhandler_m,
dofhandler_s,
face_q_collection);
unsigned int
ndofs_bi = interface_data.ndofs_biorthogonal();
ndofs_m = dofhandler_m.n_dofs();
ndofs_s = dofhandler_s.n_dofs();
n_cells_m = interface_data.n_cells_master();
n_cells_s = interface_data.n_cells_slave();
{
constraints_m.clear();
constraints_s.clear();
constraints.clear();
DoFTools::make_hanging_node_constraints (dofhandler_m,
constraints_m);
DoFTools::make_hanging_node_constraints (dofhandler_s,
constraints_s);
class DirichletBoundary1 : public Function<3>
{
public:
DirichletBoundary1 () : Function<3> (3) {};
void
vector_value (Point<3> const& p,
Vector<double>& value) const { value (0) = 0.1*std::sin((p[1]-1.0)/3.0*numbers::PI);
value (1) = 0.0;
value (2) = 0.0;};
void
vector_value_list (std::vector<Point<3>> const& points,
std::vector<Vector<double>>& values) const { for (unsigned int j=0; j < 3; ++j) {DirichletBoundary1::vector_value(points[j],values[j]);}};
};
VectorTools::interpolate_boundary_values (dofhandler_m,
2,
ZeroFunction<3> (tdim),
constraints_m,
ComponentMask (std::vector<bool> ({true,true,false})));
VectorTools::interpolate_boundary_values (dofhandler_m,
4,
ConstantFunction<3> (-0.1, tdim),
constraints_m,
ComponentMask (std::vector<bool> ({false,false,true})));
VectorTools::interpolate_boundary_values (dofhandler_s,
1,
DirichletBoundary1 (),
constraints_s,
ComponentMask (std::vector<bool> ({true,true,false})));
VectorTools::interpolate_boundary_values (dofhandler_s,
3,
ZeroFunction<3> (tdim),
constraints_s,
ComponentMask (std::vector<bool> ({false,false,true})));
BlockDynamicSparsityPattern dsp (3,3);
dsp.block(0,0).reinit (ndofs_m, ndofs_m);
dsp.block(0,1).reinit (ndofs_m, ndofs_s);
dsp.block(0,2).reinit (ndofs_m, ndofs_bi);
dsp.block(1,0).reinit (ndofs_s, ndofs_m);
dsp.block(1,1).reinit (ndofs_s, ndofs_s);
dsp.block(1,2).reinit (ndofs_s, ndofs_bi);
dsp.block(2,0).reinit (ndofs_bi, ndofs_m);
dsp.block(2,1).reinit (ndofs_bi, ndofs_s);
dsp.block(2,2).reinit (ndofs_bi, ndofs_bi);
dsp.collect_sizes();
DoFTools::make_sparsity_pattern(dofhandler_m,
dsp.block(0,0),
constraints_m,
false);
DoFTools::make_sparsity_pattern(dofhandler_s,
dsp.block(1,1),
constraints_s,
false);
// make_biorthogonal_sparsity_pattern (dsp);
deallog << "Nonzero entries: " << dsp.n_nonzero_elements() << std::endl;
timer.enter_subsection("copy_sparsity_pattern");
sparsity_pattern.copy_from(dsp);
timer.leave_subsection("copy_sparsity_pattern");
constraints_s.shift (ndofs_m);
constraints.merge (constraints_m);
constraints.merge (constraints_s);
constraints.close ();
}
system_matrix.reinit(sparsity_pattern);
deallog << "\t n_dofs dofhandler_m: "
<< ndofs_m
<< "\n\t n_dofs dofhandler_s: "
<< ndofs_s
<< "\n\t n_dofs biorthogonal basis: "
<< ndofs_bi
<< "\n\t n_dofs all: "
<< int(ndofs_bi + ndofs_m + ndofs_s)
<< "\n"
<< std::endl;
}
int main ()
{
try
{
dealii::deallog.depth_console(1);
std::ofstream out ("run/output_data/Run.log");
dealii::deallog.attach(out);
HybridMesh test;
test.run();
}
catch (std::exception &exc)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Exception on processing: " << std::endl
<< exc.what() << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
catch (...)
{
std::cerr << std::endl << std::endl
<< "----------------------------------------------------"
<< std::endl;
std::cerr << "Unknown exception!" << std::endl
<< "Aborting!" << std::endl
<< "----------------------------------------------------"
<< std::endl;
return 1;
}
return 0;
}
JobId thinkpad-lsx.site Mon Aug 15 12:35:24 2016
setup_intial_fe_distribution, wall time: 0.063519s.
DEAL::Nonzero entries: 6347893
copy_sparsity_pattern, wall time: 56.5698s.
DEAL:: n_dofs dofhandler_m: 240
n_dofs dofhandler_s: 131769
n_dofs biorthogonal basis: 1089
n_dofs all: 133098
setup_system, wall time: 60.2581s.
Setup, wall time: 60.3217s.
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 61.9s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Setup | 1 | 60.3s | 97.% |
| copy_sparsity_pattern | 1 | 56.6s | 91.% |
| setup_intial_fe_distribution | 1 | 0.0635s | 0.10% |
| setup_system | 1 | 60.3s | 97.% |
+---------------------------------+-----------+------------+------------+
JobId thinkpad-lsx.site Mon Aug 15 12:39:24 2016
setup_intial_fe_distribution, wall time: 0.0154519s.
DEAL::Nonzero entries: 6347893
copy_sparsity_pattern, wall time: 20.1535s.
DEAL:: n_dofs dofhandler_m: 240
n_dofs dofhandler_s: 131769
n_dofs biorthogonal basis: 1089
n_dofs all: 133098
setup_system, wall time: 20.7195s.
Setup, wall time: 20.7351s.
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 21.2s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Setup | 1 | 20.7s | 98.% |
| copy_sparsity_pattern | 1 | 20.2s | 95.% |
| setup_intial_fe_distribution | 1 | 0.0155s | 0.073% |
| setup_system | 1 | 20.7s | 98.% |
+---------------------------------+-----------+------------+------------+
