Dear Martin,
unfortunately I have to revise all my statements. Updating the deal.II
library to the most recent developement version did not solve the problem.
I made a stupid mistake while creating the file for the test case. The
function make_offdiagonal_sparsity_pattern () has been commented out in the
test case file. And with this line commented out I need the same time to
copy as you did. If I take into account this function the time to copy the
dynamicsparsitypattern becomes as big as before. I attached the corrected
testcase file to the post again and also the two log files one for each
run. Can you run the program on your machine with the modified file again
and figure out whether you get similiar results or not? If you observe the
same behavior we can be sure that the problem is related to my function
make_offdiagonal_sparsity_pattern
() and the way I create the sparsity pattern.
Best,
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(2);
}
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 << std::endl;
for (unsigned int block_i = 0; block_i < dsp.n_block_rows(); ++block_i)
{
for (unsigned int block_j = 0; block_j < dsp.n_block_cols(); ++block_j)
{
deallog << "Memory used by block (" << block_i << ", " << block_j << ") = " << dsp.block(block_i, block_j).memory_consumption() << std::endl;
}
}
deallog << std::endl;
deallog <<"Nonzero entries: " << dsp.n_nonzero_elements() << std::endl;
deallog << 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\t n_dofs dofhandler_s: "
<< ndofs_s
<< "\n\t\t n_dofs biorthogonal basis: "
<< ndofs_bi
<< "\n\t\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 Tue Aug 30 14:01:00 2016
setup_intial_fe_distribution, wall time: 0.581621s.
DEAL::
DEAL::Memory used by block (0, 0) = 50141
DEAL::Memory used by block (0, 1) = 5901
DEAL::Memory used by block (0, 2) = 5901
DEAL::Memory used by block (1, 0) = 20909301
DEAL::Memory used by block (1, 1) = 344734525
DEAL::Memory used by block (1, 2) = 20909301
DEAL::Memory used by block (2, 0) = 86901
DEAL::Memory used by block (2, 1) = 86901
DEAL::Memory used by block (2, 2) = 86901
DEAL::
DEAL::Nonzero entries: 53097165
DEAL::
copy_sparsity_pattern, wall time: 0.386629s.
DEAL:: n_dofs dofhandler_m: 240
n_dofs dofhandler_s: 871215
n_dofs biorthogonal basis: 3615
n_dofs all: 875070
setup_system, wall time: 29.6744s.
Setup, wall time: 30.257s.
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 39.9s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Setup | 1 | 30.3s | 76.% |
| copy_sparsity_pattern | 1 | 0.387s | 0.97% |
| setup_intial_fe_distribution | 1 | 0.582s | 1.5% |
| setup_system | 1 | 29.7s | 74.% |
+---------------------------------+-----------+------------+------------+
JobId thinkpad-lsx.site Tue Aug 30 14:02:48 2016
setup_intial_fe_distribution, wall time: 0.657087s.
DEAL::
DEAL::Memory used by block (0, 0) = 50141
DEAL::Memory used by block (0, 1) = 5901
DEAL::Memory used by block (0, 2) = 180501
DEAL::Memory used by block (1, 0) = 20909301
DEAL::Memory used by block (1, 1) = 344734525
DEAL::Memory used by block (1, 2) = 21246729
DEAL::Memory used by block (2, 0) = 261501
DEAL::Memory used by block (2, 1) = 424329
DEAL::Memory used by block (2, 2) = 86901
DEAL::
DEAL::Nonzero entries: 53353179
DEAL::
copy_sparsity_pattern, wall time: 15.5338s.
DEAL:: n_dofs dofhandler_m: 240
n_dofs dofhandler_s: 871215
n_dofs biorthogonal basis: 3615
n_dofs all: 875070
setup_system, wall time: 44.9071s.
Setup, wall time: 45.5652s.
+---------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 55.7s | |
| | | |
| Section | no. calls | wall time | % of total |
+---------------------------------+-----------+------------+------------+
| Setup | 1 | 45.6s | 82.% |
| copy_sparsity_pattern | 1 | 15.5s | 28.% |
| setup_intial_fe_distribution | 1 | 0.657s | 1.2% |
| setup_system | 1 | 44.9s | 81.% |
+---------------------------------+-----------+------------+------------+
