From fc4d48b4ca8a7bd66098c0ef8c589dbb4c4eb5b8 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Fri, 25 Jan 2019 22:18:17 +0100 Subject: [PATCH] Add internal::MatrixFreeFunctions::DoFInfo::get_dof_indices_on_cell_batch() --- .../images/dofinfo_get_dof_indices.png | Bin 0 -> 6257 bytes doc/news/changes/minor/20190130DenisDavydov | 4 + include/deal.II/matrix_free/dof_info.h | 24 ++ .../deal.II/matrix_free/dof_info.templates.h | 67 ++++++ tests/matrix_free/dof_info_01.cc | 203 +++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=1.output | 34 +++ tests/matrix_free/dof_info_02.cc | 205 ++++++++++++++++++ ...h_mpi=true.with_p4est=true.mpirun=1.output | 34 +++ 8 files changed, 571 insertions(+) create mode 100644 doc/doxygen/images/dofinfo_get_dof_indices.png create mode 100644 doc/news/changes/minor/20190130DenisDavydov create mode 100644 tests/matrix_free/dof_info_01.cc create mode 100644 tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output create mode 100644 tests/matrix_free/dof_info_02.cc create mode 100644 tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output diff --git a/doc/doxygen/images/dofinfo_get_dof_indices.png b/doc/doxygen/images/dofinfo_get_dof_indices.png new file mode 100644 index 0000000000000000000000000000000000000000..8366cacc63304bf9d7684113dca14ee7c87603d7 GIT binary patch literal 6257 zcmb_?nKzpmf!z3aXAkGs~r&pQ2_*!#2hK6{_=M>?wHqzt400FbMzDd__M2!YEf#CW(7 zy7RRpTt(udX7U;UD7yb-P&^L>BLI->%LGS6Z5ft(~Y90C|jFCkEQJQ`;Xfg<>hU z6q3zB@E}n8o0x~1#PFK{#4Ac6AAqQU;OR-620#`*@N~e&dIfmI2Ryy!w>}EM)6bH< zL4a8tD=jE11_0lHg((5gq=3?)#}Uea@dMzdoo2TLFvSNvP>0&70iWxDb~M#>BtT3C zJkSdZ<^u4&fu|^Tb`KyZ8Mvu(0G0m9Q%<$PjYBFKDc!&!q3mx)%=d~23gy3ZAFY0i z=cbs|HLDZ|f2U_My^y~M<@dv{0FV_!i{tk2(gRIVjz)_`)lyjUZM6|%*{rQsuGWT1 zTp$3j^u~AiidUeD0WL)de{q#_7lZ%I;(C_bX`~&J5>f+XZO%XsUH)LBoQ-Iho?csB zozm!5Fth51`dry}TDC(kpa1fe{(XA*$xU+9UyD*oB`&aX#Wc}zneNrJ5B)}vK_ z*Imc{W*%khR<~|2=DxY2NAoI5Ddqj4P##;j!bh&2yCYcJ( z!~md_<q5cbX;e=k<3f3EMV!%bzuBJzUmit#!-{n#3}hW)%pOI2D>?X#wL=Y` zCXxgECG4gin`VrtnxKA!A=pTb`9Q|?nnb9%N@I)^3J?nT9n!2S80U}FlNlqcQEVs_ zX9!Zsx&F=ejbKVNT~_Kh*J=it_hQ)_-%6fQ&P2;9_C09-ahpiE`}qU@-Y$j8hrA@Q zHUpJxrDP?&44)sD;x&aI+YD8@cIsTG2mI7FytNW@9- zg7iA?vBef@R&zE||LVSNEhHFftj5beK@)vLfXX{ev0I#jCR|bIUIXHqdS$NR2+fGf z$V+|pP2QqZ5j6((*e^4-PDmn&I9@OhDtI-0wQE&)m2s8zz>GBA8uBXl%wYMxe!rsz z{TkJp$QmImNl4!4Q;Pn4VYVT)Xtb$9S$3^{g}5kowSimne1$2%w5a346ro02Y>n<8b3i?SsYg*AksRr-f}!xyQB4#h5Rm)gTBncPdUuOfW7{E72=a-EI*p-z`YbZ&7YCx$viu@6)(veoxgt z8cmSCVK8klUqD(wr!6CST+UzYl*gFeq~zZAvH+8m(P{zD&~h$V@NC&QXYvuaG`xIq zUT~~RW=vQ`#z8nshQ2M`Bl^l&npvQ4_NSiWlP-v2o>jb6;!;pd!?5(Hag`4$l9Ui7}`{lrPg$;zpzy`_yw(bsOs=gODTlG5rM{5qI3IFT%xDkC^4w2``zI+?my z_tezGl()99wjQRATs6%%!PGXt3^tW8$ul-GeNv@TmR@#uqN%K@fV!YTCswB#;>nZQT#J0rxH2QqFlp!Jyou(*ZHf8pGUF~{P z^0DMX3OWrV=hJhPuqj|~JVZ4#CA?C5V}PsMA>AN#RMq~LeZ%X(!?4J3&42;^L|Oh+ z#@=~U%c{nk<5y3Xyhe`GQ>#)tcME4LFk%=7Oz2$ZSCLGv3=`qk-?~tW-}Tmp<*sMO zKX6(s_0?6>^|!CtU>C3vh=G<(ey+=NFFv)}yCr)f>zcVDyCN;VAhcj&`Bl8V4pldx z+vS{L8N8_I$BwnddSQow0{C5!JzluKyuaDCi!-lxtZh+plC!WAY0eTKVXJ zvVgUYGS+9P9Yh)H=|&Pd2JJu*{K)sm#5W?Vb*soOU$T3_p~27)xhR(~V}@bDL{azQ zImOqDO6N)mN(-s>g1QnFd>W$Mvc>$4(lg?XB8_sr{B|Z^$BomC)t_rJN$_)*i@ki_ zy7+tdFmyS6`TC0Ytt0`|6Bo~HbmviU0bLP44X-0wz_ zdUi5*+Or7^iw%cODVlTN(K%4#*7o|i6Voi}$QdltZs^?{7oxq)VqAyNWh*AO*KGIr+{2zv< z_b0iDkoGXBUE5%IzeGYw2G8K`*Q)+mgm~IsCU%#mlV<5-ZC;2bl#khXy`|Jq0xA8< zbCz_4mVhVZF=456eMBwyi*enOU1(ytZwsfFwv)d_M}I&^4ajNxx)JGGu+R3~@G^4tDC1Ji>G(U zC1>8g_by!`8fAIQY|kvqqZZGfbl8dU3<;U} zmP0S+)1K4gMY?%k*L#=ffu?Eqdzi9j-Db-p$gcdf;(X}1oY~df)ZUz*W4n9S;qaX5 zDPeHn&{fQFG^HXzbaM1S2zkhM#zR>NiDfz0%gx^fibJ+ToEf*Sa;^r-Z_x?(o{eu$ zqKA`)vl9;zW7}T5%e&la-=?pa?d|Nn)B43jSkU6=@x^Y3`P-M}djJ4iPXRy@ z7n9y{{%F~w0swfxBP~N^EEYEf$jZw4Q?=t}W7`2NfWfpw02B(4cUZ!{kpcc zhC52H@V{!ao+Xw41_Obg}wSyhl zj=}cN@34?8EV#9`6`PrfWnsZCE@CY$v3`Eowl=Js919D};^JakTbre&rJS5xW@e^2 z5aU0P#mC1dBqSszCMG2%;b2cmNl8sjO-oBlPfyRt$iNYXlWcZ&c1}*t zr%#`9b93|Z^78ZZaa0u+78Vs16&Dwml$4Z~mX?*3m6w-SR8)Na{JFBS5{X3OBv@Tt zT~kw2TU%RKS65$O-_X#|*x1U0q$>-Q7JsJ-xlX zeSLi>6za>DFa7=f0|Ns%H5eQmM5EC|Lqp%beH$Jg9vK-K9UUDT8yg=VpO~1KoSdAR znwp-To|&1Mot>SVo134XUszbcf&BgZ_obz!<>lp-m6g@iRh&Aludi=xZ2b7~V{>zJ zYinzJdwXYRXLolOgTd_W?d|XH|NQy$;Nal!@bKv9==k{f*S&C!*9zR#BJHzY|C)k@W3o|E+dxCG?&kxW6HvnnT0*+qEp*1Es|S9RZgQ)t z-*lPy=Nnbtuxf@H2Gi9z+DT@7c^pZ9(s4|Pu_*>M>igWeAaLjlB}<~zK{dO4d8Qcn z5YHinQPPgaXk{J<|OFVVBH~kI3`?+@CU^QcEHoXctqoRw;9h7z_vg z_XR2kvLKdPmt}EXtHlCwnmM~`$2j8TUoPfm^)L@V&~Eshm0H#Cha z2=XliQL|3}G`LU$L%VNpyk;cIw>TvB{9u;airj(O5dTk!Lwx`{UPNzFR&bise9%mX z66uRf?+acoMRIOSRh3h!Kas40Gj4W(U+6K)=F9z-eSy|18aRh^7VySwID(EifBq{u+LsnvTu&jNBty>7)nucHe)@WPs9a36upn5D_Tf$%Rl{<{{!&A(8+V#Zs@Zv1gX{)G40j+XRt zKhm)Yf+$T95aUYuI;OSj)!8b%eKleIWKZgyUR8OJxyfSkp2#|Ze&fKCc{Q`K^T!fUue#ZLeU-M(xJt|ksr3d-a!-$|HZv#WV8`hulPXV=oeYO4I4pNsOsTfd z9}mjF)Kf)N{dr&Z&QR`XZpgsv3wTG_8r-%B&^&_}^8^ytO*&>sGo_os4K7j2W3)zh z^9Zj?>v>xD@fPtBTx+d%)2TTLEnN+~&Oz0a2`$gbCMylFf|hGf+g^k-sOJMb4KD8&4riiM`x@_)^A!DZlkCgK)%{NNhABgj%SA>pJ2o)nk)G zQ<%}7BEt#m;mqx8|7N}iGVCgKK-W&K>Z4__zBE5sRC#+5s+Tqa{QI%Io;?kRiWDR8l@Ot?ZQ8#nt zIEGQS6n#nWQ5OH5*Qb173q;#kIkZsXCEqd=D@%jaW0RMNMXtO!vkm0wY5BEi%3OYN zFnjrXxksS?1PjU-5nDmyV~J>V|MH|v4hdbAz2qJiCoC$i{qA^4q7VJDv|pkGm%d}` zjfw7Er0lk=>@!;2_#`>hFe7g1ZHP90kC>u%JuqT}+#%DKvaIYUwScAj#)tV3Gb;Z7 zYv5mckMQ=N%jPEed$m>(6qUFwJCUQF6oY=5;kgz#LKY_A4`6=@q`8jdx%b-KWaas- zJj+M0c#>th$*JXD>PFYnDeOf)pKt5 zE|8Ff8cEUgPY=Yr^bWotAKOAvlZ1w}s2ZBu=*1lhq`MlA&|K zpbBNbaW6+iLp133qB)Dc*qVxx`p zb&V)^G_bjae11pMV}yYUl7vh`dK6aez@{(Ms6?G%=|dAO_BKgIW}8D|gsYl}oG;Ds zZXpxOgN?4tt5ExiOJ;p`yFr|nIea`pZ@TfWmg!qtektz0DrLq;4RfrPlnPg4a1P1M zKk15uOC9-^3b!Y-WvWOMq;y*t%noU8sl}t-h=ZBt+V?XED^%MB#i7wgSyI=gCn=V* zWgVxt{X*1q%3*tt9Kv9GWsF^B=l&m9|1$U2Tckzz-|B{v3SD|{WrB?whOklN;&{o; zAG<(m+L*&O8v>Lfj;)+x5P7hlO56@8h#F@mktQ%uXRnJd*k8F6^p8(75kyV#IIWL3 z{TCl}-~)N-$;W`$lz++11mK_(kDmZ+y%4DpL%2A?fX;x9(k-Ed*PDeThf|qH!)Z0n zJpLgAwezwdlz}j$fc=NR#xn`%mf|@ z^SFWMDyguAv$A!JN_Xm+l(xwo#IkQ;WMO#G>3(qN`F*O*k_=WKXr&lz+Kim0hQ-F_ z;h+xlzwDL@@|fTcbT(#C<)>RHm$Z|kc&6Uk@+M^cST=Z4mk~i9pk+q>?D=*vu@*CZ z)VbQ`_tCpd6OAX2V_{KLAX_VNDRRB0pN%*G{O@TwwFgTC$OX=SAKIv?#J#Qnb!8o; Ja)qY>{{fV@wLbs= literal 0 HcmV?d00001 diff --git a/doc/news/changes/minor/20190130DenisDavydov b/doc/news/changes/minor/20190130DenisDavydov new file mode 100644 index 0000000000..4cf4ec962c --- /dev/null +++ b/doc/news/changes/minor/20190130DenisDavydov @@ -0,0 +1,4 @@ +New: Add internal::MatrixFreeFunctions::DoFInfo::get_dof_indices_on_cell_batch() to +return locally owned DoFs used by matrix-free framework on a given cell. +
+(Denis Davydov, 2019/01/30) diff --git a/include/deal.II/matrix_free/dof_info.h b/include/deal.II/matrix_free/dof_info.h index 6839cd1d35..8f80cf2334 100644 --- a/include/deal.II/matrix_free/dof_info.h +++ b/include/deal.II/matrix_free/dof_info.h @@ -102,6 +102,30 @@ namespace internal fe_index_from_degree(const unsigned int first_selected_component, const unsigned int fe_degree) const; + /** + * Populate the vector @p locall_indices with locally owned degrees of freedom + * stored on the cell block @p cell. + * If @p with_constraints is `true`, then the returned vector will contain indices + * required to resolve constraints. + * + * The image below illustrates the output of this function for cell blocks + * zero and one with zero Dirichlet boundary conditions at the bottom of + * the domain. Note that due to the presence of constraints, the DoFs + * returned by this function for the case `with_constraints = true` are + * not a simple union + * of per cell DoFs on the cell block @p cell. + * + * @image html dofinfo_get_dof_indices.png + * + * @note The returned indices may contain duplicates. The unique set can be + * obtain using `std::sort()` followed by `std::unique()` and + * `std::vector::erase()`. + */ + void + get_dof_indices_on_cell_batch(std::vector &locall_indices, + const unsigned int cell, + const bool with_constraints = true) const; + /** * This internal method takes the local indices on a cell and fills them * into this class. It resolves the constraints and distributes the diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 153e6a8236..a2a7c6ae68 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -183,6 +183,73 @@ namespace internal + void + DoFInfo::get_dof_indices_on_cell_batch(std::vector &my_rows, + const unsigned int cell, + const bool apply_constraints) const + { + const unsigned int n_fe_components = start_components.back(); + const unsigned int fe_index = + dofs_per_cell.size() == 1 ? 0 : cell_active_fe_index[cell]; + const unsigned int dofs_this_cell = dofs_per_cell[fe_index]; + + const unsigned int n_vectorization = vectorization_length; + constexpr auto dof_access_index = dof_access_cell; + AssertIndexRange(cell, + n_vectorization_lanes_filled[dof_access_index].size()); + const unsigned int n_vectorization_actual = + n_vectorization_lanes_filled[dof_access_index][cell]; + + // we might have constraints, so the final number + // of indices is not known a priori. + // conservatively reserve the maximum without constraints + my_rows.reserve(n_vectorization * dofs_this_cell); + my_rows.resize(0); + unsigned int total_size = 0; + for (unsigned int v = 0; v < n_vectorization_actual; ++v) + { + const unsigned int ib = + (cell * n_vectorization + v) * n_fe_components; + const unsigned int ie = + (cell * n_vectorization + v + 1) * n_fe_components; + + // figure out constraints by comparing constraint_indicator row + // shift for this cell within the block as compared to the next + // one + const bool has_constraints = + row_starts[ib].second != row_starts[ib + n_fe_components].second; + + auto do_copy = [&](const unsigned int *begin, + const unsigned int *end) { + const unsigned int shift = total_size; + total_size += (end - begin); + my_rows.resize(total_size); + std::copy(begin, end, my_rows.begin() + shift); + }; + + if (!has_constraints || apply_constraints) + { + const unsigned int *begin = + dof_indices.data() + row_starts[ib].first; + const unsigned int *end = + dof_indices.data() + row_starts[ie].first; + do_copy(begin, end); + } + else + { + Assert(row_starts_plain_indices[cell * n_vectorization + v] != + numbers::invalid_unsigned_int, + ExcNotInitialized()); + const unsigned int *begin = + plain_dof_indices.data() + + row_starts_plain_indices[cell * n_vectorization + v]; + const unsigned int *end = begin + dofs_this_cell; + do_copy(begin, end); + } + } + } + + template void DoFInfo ::read_dof_indices( diff --git a/tests/matrix_free/dof_info_01.cc b/tests/matrix_free/dof_info_01.cc new file mode 100644 index 0000000000..462abcac78 --- /dev/null +++ b/tests/matrix_free/dof_info_01.cc @@ -0,0 +1,203 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + + +// Tests DoFInfo::get_dof_indices() + +#include + +#include +#include + +#include + +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + + + +template +void +test(const bool adaptive_ref = true) +{ + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + const unsigned int this_mpi_core = + dealii::Utilities::MPI::this_mpi_process(mpi_communicator); + + parallel::distributed::Triangulation tria(mpi_communicator); + GridGenerator::hyper_cube(tria, 0, 1, true); + tria.refine_global(1); + if (adaptive_ref) + { + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + if (cell->center().norm() < 0.5) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + if (cell->center()[0] < 0.2) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + else + { + tria.refine_global(1); + } + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + IndexSet owned_set = dof.locally_owned_dofs(); + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs(dof, relevant_set); + + AffineConstraints constraints(relevant_set); + DoFTools::make_hanging_node_constraints(dof, constraints); + // constrain bottom part of the boundary (lower in y direction) + VectorTools::interpolate_boundary_values(dof, + 2, + Functions::ZeroFunction(), + constraints); + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + std::shared_ptr> mf_data( + new MatrixFree()); + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.tasks_block_size = 7; + mf_data->reinit(dof, constraints, quad, data); + } + + const unsigned int n_cells = mf_data->n_macro_cells(); + const auto & dof_info = mf_data->get_dof_info(); + constexpr unsigned int n_vectorization = + VectorizedArray::n_array_elements; + + std::vector my_rows; + my_rows.reserve(fe.dofs_per_cell * n_vectorization); + + constexpr auto dof_access_index = + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; + + bool checked_interleaved = false; + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + deallog << "Cell: " << cell << std::endl; + checked_interleaved |= + (dof_info.index_storage_variants[dof_access_index][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved); + auto get_and_log = [&](const bool apply_constraints) { + dof_info.get_dof_indices_on_cell_batch(my_rows, + cell, + apply_constraints); + // sort and make unique to make visual inspection easier: + std::sort(my_rows.begin(), my_rows.end()); + my_rows.erase(std::unique(my_rows.begin(), my_rows.end()), + my_rows.end()); + for (auto el : my_rows) + deallog << " " << el; + deallog << std::endl; + }; + get_and_log(true); + get_and_log(false); + } + + // if we don't have adaptive refinement, we should have + // some cells interleaved (at least that is the intention of this + // part of the test) + Assert(checked_interleaved || adaptive_ref, ExcInternalError()); + + // output in Gnuplot + if (dim == 2) + { + std::map> support_points; + MappingQ1 mapping; + DoFTools::map_dofs_to_support_points(mapping, dof, support_points); + + const std::string prefix = + std::is_same::value ? "float_" : "double_"; + const std::string href = (adaptive_ref ? "" : "global_"); + const std::string base_filename = + prefix + href + "grid" + dealii::Utilities::int_to_string(dim) + "_" + + dealii::Utilities::int_to_string(fe_degree) + "_p" + + dealii::Utilities::int_to_string(this_mpi_core); + + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" + << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "unset grid" << std::endl + << "unset border" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels tc rgb 'red' nopoint notitle, '-' with labels point pt 4 offset 1,1 notitle" + << std::endl; + GridOut().write_gnuplot(tria, f); + f << "e" << std::endl; + + // output cell blocks: + for (unsigned int cell = 0; cell < n_cells; ++cell) + for (unsigned int c = 0; c < mf_data->n_components_filled(cell); ++c) + { + const auto dof_cell = mf_data->get_cell_iterator(cell, c); + f << dof_cell->center() << " \"" << cell << "\"\n"; + } + + f << std::flush; + f << "e" << std::endl << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, support_points); + + f << "e" << std::endl; + } +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + MPILogInitAll mpi_init_log; + + deallog.push("2d"); + test<2, 1>(); + test<2, 1, float>(); + test<2, 1>(false); + deallog.pop(); +} diff --git a/tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..226a54f496 --- /dev/null +++ b/tests/matrix_free/dof_info_01.with_mpi=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,34 @@ + +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 7 18 19 21 22 +DEAL:0:2d:: 6 7 16 17 18 19 20 21 22 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 2 7 9 10 21 22 23 24 +DEAL:0:2d:: 0 6 7 8 10 21 22 23 24 25 26 +DEAL:0:2d::Cell: 2 +DEAL:0:2d:: 2 3 7 9 10 11 12 24 +DEAL:0:2d:: 0 1 2 3 7 8 9 10 11 12 24 25 26 +DEAL:0:2d::Cell: 3 +DEAL:0:2d:: 2 3 4 5 9 11 12 14 15 +DEAL:0:2d:: 2 3 4 5 9 11 12 13 14 15 +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 2 7 9 10 18 19 21 22 23 24 +DEAL:0:2d:: 0 6 7 8 10 16 17 18 19 20 21 22 23 24 25 26 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 2 3 4 5 7 9 10 11 12 14 15 24 +DEAL:0:2d:: 0 1 2 3 4 5 7 8 9 10 11 12 13 14 15 24 25 26 +DEAL:0:2d::Testing FE_Q<2>(1) +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 2 3 5 6 7 8 +DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 5 8 10 12 13 14 +DEAL:0:2d:: 4 5 8 9 10 11 12 13 14 +DEAL:0:2d::Cell: 2 +DEAL:0:2d:: 6 7 8 15 16 17 18 19 20 +DEAL:0:2d:: 6 7 8 15 16 17 18 19 20 +DEAL:0:2d::Cell: 3 +DEAL:0:2d:: 8 13 14 17 20 21 22 23 24 +DEAL:0:2d:: 8 13 14 17 20 21 22 23 24 diff --git a/tests/matrix_free/dof_info_02.cc b/tests/matrix_free/dof_info_02.cc new file mode 100644 index 0000000000..a980641a68 --- /dev/null +++ b/tests/matrix_free/dof_info_02.cc @@ -0,0 +1,205 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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. +// +// --------------------------------------------------------------------- + + + +// Tests DoFInfo::get_dof_indices() for multi-component system + +#include + +#include +#include + +#include +#include + +#include +#include + +#include + +#include + +#include + +#include + +#include "../tests.h" + + + +template +void +test(const bool adaptive_ref = true) +{ + MPI_Comm mpi_communicator(MPI_COMM_WORLD); + const unsigned int this_mpi_core = + dealii::Utilities::MPI::this_mpi_process(mpi_communicator); + + parallel::distributed::Triangulation tria(mpi_communicator); + GridGenerator::hyper_cube(tria, 0, 1, true); + tria.refine_global(1); + if (adaptive_ref) + { + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + if (cell->center().norm() < 0.5) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + for (auto cell : tria.active_cell_iterators()) + if (cell->is_locally_owned()) + if (cell->center()[0] < 0.2) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + else + { + tria.refine_global(1); + } + + FE_Q feq(fe_degree); + FESystem fe(feq, dim); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + IndexSet owned_set = dof.locally_owned_dofs(); + IndexSet relevant_set; + DoFTools::extract_locally_relevant_dofs(dof, relevant_set); + + AffineConstraints constraints(relevant_set); + DoFTools::make_hanging_node_constraints(dof, constraints); + // constrain bottom part of the boundary (lower in y direction) + VectorTools::interpolate_boundary_values(dof, + 2, + Functions::ZeroFunction(2), + constraints); + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + std::shared_ptr> mf_data( + new MatrixFree()); + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.tasks_block_size = 7; + mf_data->reinit(dof, constraints, quad, data); + } + + const unsigned int n_cells = mf_data->n_macro_cells(); + const auto & dof_info = mf_data->get_dof_info(); + constexpr unsigned int n_vectorization = + VectorizedArray::n_array_elements; + + std::vector my_rows; + my_rows.reserve(fe.dofs_per_cell * n_vectorization); + + constexpr auto dof_access_index = + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; + + bool checked_interleaved = false; + for (unsigned int cell = 0; cell < n_cells; ++cell) + { + deallog << "Cell: " << cell << std::endl; + checked_interleaved |= + (dof_info.index_storage_variants[dof_access_index][cell] == + internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: + interleaved); + auto get_and_log = [&](const bool apply_constraints) { + dof_info.get_dof_indices_on_cell_batch(my_rows, + cell, + apply_constraints); + // sort and make unique to make visual inspection easier: + std::sort(my_rows.begin(), my_rows.end()); + my_rows.erase(std::unique(my_rows.begin(), my_rows.end()), + my_rows.end()); + for (auto el : my_rows) + deallog << " " << el; + deallog << std::endl; + }; + get_and_log(true); + get_and_log(false); + } + + // if we don't have adaptive refinement, we should have + // some cells interleaved (at least that is the intention of this + // part of the test) + Assert(checked_interleaved || adaptive_ref, ExcInternalError()); + + // output in Gnuplot + if (dim == 2) + { + std::map> support_points; + MappingQ1 mapping; + DoFTools::map_dofs_to_support_points(mapping, dof, support_points); + + const std::string prefix = + std::is_same::value ? "float_" : "double_"; + const std::string href = (adaptive_ref ? "" : "global_"); + const std::string base_filename = + prefix + href + "grid" + dealii::Utilities::int_to_string(dim) + "_" + + dealii::Utilities::int_to_string(fe_degree) + "_p" + + dealii::Utilities::int_to_string(this_mpi_core); + + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" + << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "unset grid" << std::endl + << "unset border" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels tc rgb 'red' nopoint notitle, '-' with labels point pt 4 offset 1,1 notitle" + << std::endl; + GridOut().write_gnuplot(tria, f); + f << "e" << std::endl; + + // output cell blocks: + for (unsigned int cell = 0; cell < n_cells; ++cell) + for (unsigned int c = 0; c < mf_data->n_components_filled(cell); ++c) + { + const auto dof_cell = mf_data->get_cell_iterator(cell, c); + f << dof_cell->center() << " \"" << cell << "\"\n"; + } + + f << std::flush; + f << "e" << std::endl << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, support_points); + + f << "e" << std::endl; + } +} + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization( + argc, argv, testing_max_num_threads()); + + MPILogInitAll mpi_init_log; + + deallog.push("2d"); + test<2, 1>(); + test<2, 1, float>(); + test<2, 1>(false); + deallog.pop(); +} diff --git a/tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output b/tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..1da9dc63ea --- /dev/null +++ b/tests/matrix_free/dof_info_02.with_mpi=true.with_p4est=true.mpirun=1.output @@ -0,0 +1,34 @@ + +DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2] +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 14 15 36 37 38 39 42 43 44 45 +DEAL:0:2d:: 12 13 14 15 32 33 34 35 36 37 38 39 40 41 42 43 44 45 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 4 5 14 15 18 19 20 21 42 43 44 45 46 47 48 49 +DEAL:0:2d:: 0 1 12 13 14 15 16 17 20 21 42 43 44 45 46 47 48 49 50 51 52 53 +DEAL:0:2d::Cell: 2 +DEAL:0:2d:: 4 5 6 7 14 15 18 19 20 21 22 23 24 25 48 49 +DEAL:0:2d:: 0 1 2 3 4 5 6 7 14 15 16 17 18 19 20 21 22 23 24 25 48 49 50 51 52 53 +DEAL:0:2d::Cell: 3 +DEAL:0:2d:: 4 5 6 7 8 9 10 11 18 19 22 23 24 25 28 29 30 31 +DEAL:0:2d:: 4 5 6 7 8 9 10 11 18 19 22 23 24 25 26 27 28 29 30 31 +DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2] +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 4 5 14 15 18 19 20 21 36 37 38 39 42 43 44 45 46 47 48 49 +DEAL:0:2d:: 0 1 12 13 14 15 16 17 20 21 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 4 5 6 7 8 9 10 11 14 15 18 19 20 21 22 23 24 25 28 29 30 31 48 49 +DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 48 49 50 51 52 53 +DEAL:0:2d::Testing FESystem<2>[FE_Q<2>(1)^2] +DEAL:0:2d::Cell: 0 +DEAL:0:2d:: 4 5 6 7 10 11 12 13 14 15 16 17 +DEAL:0:2d:: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 +DEAL:0:2d::Cell: 1 +DEAL:0:2d:: 10 11 16 17 20 21 24 25 26 27 28 29 +DEAL:0:2d:: 8 9 10 11 16 17 18 19 20 21 22 23 24 25 26 27 28 29 +DEAL:0:2d::Cell: 2 +DEAL:0:2d:: 12 13 14 15 16 17 30 31 32 33 34 35 36 37 38 39 40 41 +DEAL:0:2d:: 12 13 14 15 16 17 30 31 32 33 34 35 36 37 38 39 40 41 +DEAL:0:2d::Cell: 3 +DEAL:0:2d:: 16 17 26 27 28 29 34 35 40 41 42 43 44 45 46 47 48 49 +DEAL:0:2d:: 16 17 26 27 28 29 34 35 40 41 42 43 44 45 46 47 48 49 -- 2.39.5