From c0c9940bc2916998fd2cbbd47bd9d42a09df1f06 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Thu, 23 Mar 2017 17:48:47 +0100 Subject: [PATCH] add a function to get DoFs that have support only within the given domain --- .../images/extract_dofs_with_support_on.png | Bin 0 -> 6912 bytes doc/news/changes/minor/20170324DenisDavydov | 4 + include/deal.II/dofs/dof_tools.h | 33 ++ source/dofs/dof_tools.cc | 104 ++++++ source/dofs/dof_tools.inst.in | 7 + ...f_tools_extract_dofs_with_support_on_01.cc | 137 ++++++++ ...ols_extract_dofs_with_support_on_01.output | 4 + ...f_tools_extract_dofs_with_support_on_02.cc | 185 ++++++++++ ...ct_dofs_with_support_on_02.mpirun=1.output | 10 + ...ct_dofs_with_support_on_02.mpirun=3.output | 20 ++ ...ct_dofs_with_support_on_02.mpirun=5.output | 30 ++ ...f_tools_extract_dofs_with_support_on_03.cc | 143 ++++++++ ...ols_extract_dofs_with_support_on_03.output | 3 + ...f_tools_extract_dofs_with_support_on_04.cc | 327 ++++++++++++++++++ ...ct_dofs_with_support_on_04.mpirun=2.output | 1 + ...ols_extract_dofs_with_support_on_04.output | 1 + 16 files changed, 1009 insertions(+) create mode 100644 doc/doxygen/images/extract_dofs_with_support_on.png create mode 100644 doc/news/changes/minor/20170324DenisDavydov create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_01.output create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_03.output create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output create mode 100644 tests/dofs/dof_tools_extract_dofs_with_support_on_04.output diff --git a/doc/doxygen/images/extract_dofs_with_support_on.png b/doc/doxygen/images/extract_dofs_with_support_on.png new file mode 100644 index 0000000000000000000000000000000000000000..9d6c3328e90a18a25189bfc024329968a88445a0 GIT binary patch literal 6912 zcma)91yodDx4#GiBBcV-Atgvi2+}bs4brKA0s{ilT?5Dff^>IEw+uO?z#uS4#}Lv2 zG9Wcd%^Sb}``%aUz4g|*C+*Vg_X6xk4ssaMBI=edA*g04O zfcJd1E?7@@n??q+bgZZr0ZCGG(xE11)lrOsP{nicvfZFmjbO{4qt@=e_4FwL10pX1 zKPCneO|8vG`krK#Xq~MfCbTf({g?9jji1s50 zX_+GPt#6?NpZizV1>}QCIGh1W;sz@=w@Y?B;8&ou^n-g{L>&O0_ry(Npi3jO>mg$} zmSRUS-3$-%4zJ5S{;?V{gbI-NjZypv$Unn_WT$cI0(k_0#n9_lYru1Uz=Atya{_>5 z|4R4615A_HsPXdR0oFU#5l;YHDWG!H;Qdp;=pjI5qxM+>nB@l^su9Q^+?0dE4`bN6vq+d72;B zGB>xezCNq^S<%$;i($YOxW~N9@Y3OTp!5aivcq&5_??Xpw znRqP(C*_$$yR>bgiO%E#J^Ib$GR_odmc%*T_uGSS5GWvs?^v?UN8r7BZ^Rx$d^bJp z7hAV7!R;sx*1-rW?R#qRkCcUU-s`dID>EI+xL%V8H+$9^FNFYvLoT50NZ4KQgQls&kAmSwQ?8Ah^5KP{ObCNUM5;Bf9q?7EyaASENI|i*LPYXkmKvnNPIn1cdwGPqMyFjpc20=^5pesofpCF;HOKhwSv778~~=3}3}9eYQR(k}w^S)Ai`BuJRMCGMI^ zU6J1S?eS;hZ*|zW`N}dymFd|N2j{Kq>xm?i_*fqx!q&6bd)GzQ8P?ejO-Zs}$-5N& z(p|l;^W}{y-3H}`=mw#6ny`X?VTR6PX};dAN3kY~RrxTT8tw4XJE9&L5+k8BDr&`5 z1<4ciws&nYb~WGKTW>8UGi#=PH5?z>R{Ea$9sk&ilnhE8IuHTAPnStcL<^?#qOD0& z|426xMK_@LP+Vv?Z7ZE3-AUhw7sl&P9sh+SSv8p}nSobOudMuI`BM2!{UiM$z1Xss zx&hVC^hI>n%5F@iROpvyKi}15&;^%Cm4=stOP%!OiUN%QU)>$Ef>&toU$zlHX5ZIofTPw-5KKfY?rN@IiUo;3vO`>J&K5mR0|mr zNR<`HWawW+z}HpXPh4KC_>P}sXEtQ^?3JQw&|+vibofHm zy>;kVzSuhoU(H^_a~xJeQzmemPclVbJwo$w(t!dcZz)YQ=iwbM z&KMG(9;P1fJ;5=tv4~la885Tup)#+AZ^DoG_D4P(eImtt8!hj+%i`MA$z2vT)rCi!7sbNgmgFJ+3()eC)o%K*j<>} zd!c0{XuRBLc?Qj9kfn}%UrmL&n@Sxtdw;=ugCX4MOK18*(x_kSIKE%j6W(xQw1o?DT_OCGr85j;ri_FMcvHGW`!({j2QGPkh$5F*lYm z%w&@Dp3#1lD?}44o=jjS!8eK-hZ)Z@bu-RQUzxUQ4a<&A?=1OVD4|#=vP7kgeN>P1 z2JPmzLnaczhojNi$=Lyc0>2k5Pz5l*t=GuqHxy3$t(5L7<4>_sTT<;(NVCIc_rrSB z#vQ|E!;Bh9bH^R`9qFsY39W&)thS8m7Bl-6v@g8T&jz#2kRy&@*y-qBe7V!LOIL&H@9Agm81xnqdUpjI{aZ`Z2yu3e67j8DT3%~+sbeBAUKmZDME7&ep03zQ7!9KKr zL;-+21i(VDLD*(p}MMOk2G&F2%Y=VM<($mu$8XCH~x`u~`H#RnK zTO0@r$D!k@W;YM?558fs=*18$th)w{#?An^-2f|G9ds}1o=$vt)0KiB% zau@g`A{Ic7VW|?aLC+5$5bRUn4C1<-;(M?e4R5k^|DJ?B6D=RB6FR!SmsI084s;a84uCA%6sjaQ8tE;Q8ug7Wd z)2C03jg2rEtf{G~xw*NerKPpCwXLnKy}ca{hyMi*PRrfh-90@$y}iAkKY#A)>+A3D z9~c-wAP|FtgI~UU85$bGnZd}&2oi}L9UcAp_3PN!*!cMP#KgqpCbF4CdFbU%!9^?xVXH$yt=x=S=P5Alq~N4#-j2>LC<^s$Kng$!pD?hE0VFP-GN5o?`$PS&ds~b z09c}WC2t5LNcP2Xa%;h~{bk%hmYCF=Hd?Yc7CV&`&JAee6u@D3SvX(M(J~Ov$s++? zH<4^ps}QDC`RHIwCe1N7`^cCrAw_P-@yE2>0vs-Kd&i#%6evxTKCG^hlIxIhY8;R{ zS`07D;o@l<5Yo^&MJksbN!|t3u@GoCjqaQiMo131#=hQ5i)wR^Nk$)Ph0JO;cH-3+ zV;*XTHo2wrBjn%X-QLSqCBJ0Q{fR-oW-1!X?*Q-ZI~YrRs<9f%NbvF)E)Fn+`jEaZ z&unZxMH`(nuo~--FKRSorLH#~+oHz}2&F)4Dom90tdOq&-aA5K))L$E%a?068U-2D z#DBnpwmmz0=Q~fW?YE(Dj9`@wbK6XZgTP`f=Bn5u&2#KV#=u(KgCeJ9#WT+bPjZ-c zXa4&4Ybh;6Ct?w*+=&nqhucaV4^ol(F$!VgV zn6y~4JKETWeDuh=wo6>~kxEOEQi~YTb>mka?NWVNgnHND;j3iyyusdF9x2Hj#Uji* zXyRh)^HcN#c)FcxJT0|&M7kY)<%1G}NcfL+^kPj(O}j25zR$vWsi}eVGWRrU{yJIq z;q#!Vt@CETPJSi^0 zY0=bg|78Ys56KGYN=hGp)E4ExE&RtX{Xvc<|Ev zZ47QEKOLyWuQ_7!K#nJ+Mb|saG-Y`B%B)whRr93A-7|Yrne{H>HoG?6YhLYj5Gudz0jW$XI!iU?Cp(6j0}Vfhr~ul!l$5N54I`f-+7Py<$f!t67)EHjUaw}qMTC!lhg1`~-!{e92Q2#LnH~Oz@^8eshNJAUnS&vf5gWmnqCqr1uhW#KtdBd=EKdn&GAS{|x{LR-e?K+(&y=9wiU zsJ18>0_C^2B{R8K8SX3Edq0Das^xT(7NJrDm!`e;FTScpup)C}Xv)T|i7)v*E*a7ibiHj{7gco6MZCAXseL@gzW z1#WnLi-S_8A=Tx!TL=LZ z9<-KKdd^oOysEvgq4m8`H6x<>FUrkJUD~JNF-Px0D*pLsu8;@AaZvweZh&j$Xbs zHQT+9v7>6SGeczo-9w1~_hhCA9CsU5)!Lv3S=;x!wvf{cK4RrloB!e4% zQQtvMKDYUowqBa0e&Q|Wk9|w&bv!I-rA4X4Y#6tj#xgo$G9?}wga2x&zF)w-RI0Q0 z^5^p)(ggmY{&m90C!>>c6lx3;)K)T+6>jWU(lhcCz7I>6UlF?}SmiuqEh7t+u zk7}1OOMO3{%D&MuWjhUhP_#!k&T>Rrmn#LLU#>rUy;T3>(~8LX%kR0TjN`SJ2cc~L z2jlx+a}kx#Azr;t_Rn#_+&w;Nv0Ktu2f2iaU!Go)(7}aR)USLDk9@uL?*_?wYad<0 zSB>JNMGSJsvV$qfs_P8#W$KAXuNMR(?<{?b(NlFL{Bmql=C4Zkm^XFya6wrTnmF_v&Nxz+#@d0x-B~@4F zB-!p1kG#q|O54j`P|+W7rN10zeG*LkH;*CzSE~4X`1ohvfi=D{p<~s$-qJ4k7f^_s z5RP<=$2zB&?Sar$Wczkjd3Mc&j4hrq)Q#-=JlW?poDai5#H|E-+kjd+8bah)uTKW$fl1~f~Qpxwpb%u0I z8KZNwQ-~^u~dDD z<9ya~kWgIlzx9Zq& z$d%fRCj98Rc_Gv_@5okwa?Y9wxyQ}8F`u~`_1&2Q8peGxbBzw~kLgPj2wj0v@8U~X z;Qyg1iCc)eCEZR3Ubx+t$x?R@7)A2B^PGWiXoN{}_2_AQy`}(zno}CgpqQzlN_L$7 zwA!m<0dqB%*!05A0AXK6iHw_(1T|a3=l|bfTrAjs5%8&^zJ`c}(6T}<0K`Sg+rw@% zc0d1ayaZR*4$CZdTxi&~59uz$S$j*TcDaLbbwsfJ{^dACmIB4+1zxDXrVZG9G|3q` zUka|V;R@0j8F#r4mIk>-W^uMm9y*}K&&-C~Rg`o*L-~r^{AK@WLk=;HlCbZ6Og7WG zS|5`PbL{~I+*om8-cfTpD#3e!{Zf67%1s#Ii{&OziQnC)pIbxKUvafvC*Jcb#o{h} zuO+t{o3uMd7f_ki)NQ9o{TXLDW>TMjqkA#rDGpF<>#>gkB0qt#Lm;4le^=R<+dTCO z4p#16TILxcu?%yDUcCPkFh<#=!~fTe@WPWa1pnpDceZ5zlTrQA7=?0z&+-e+e6Li$ z1`apRn$#;`mAOAf3|3aPk$YA znxyn%Ts*MiY#9ez+4#<4*p++9nn*R8_$ZXOUgWo={chq*vFHH%Y#YD31sP^&? zU0+puZc-TFY7wR$>d|L>X1NRwS%;dY+_V~*D$8!?B(ix$=zqsqFJO zgoe@_rN%AvmhisJ6 z-U~V_G#Yo7&fG5ExWm+8JtyR3{2^B)XF=c@{+~==R36C3yGBe#DR2At +(Denis Davydov, 2017/03/24) diff --git a/include/deal.II/dofs/dof_tools.h b/include/deal.II/dofs/dof_tools.h index aa4d0ababb..d0eead0241 100644 --- a/include/deal.II/dofs/dof_tools.h +++ b/include/deal.II/dofs/dof_tools.h @@ -1436,6 +1436,39 @@ namespace DoFTools std::vector &selected_dofs, const std::set &boundary_ids = std::set()); + /** + * Extract all degrees of freedom (DoFs) that have support only within cells, + * for which @p predicate is true. The result is returned as + * an IndexSet. + * + * Consider the following FE space where predicate returns true + * for all cells on the left side: + * + * @image html extract_dofs_with_support_on.png + * + * This functions will return the union of all DoFs on those cells minus + * DoF 11, 13, 2 and 0; the result will be [9,10], 12, [14,38]. + * On the image above the returned DoFs are separated from the rest by the + * red line + * + * Essentially, the question this functions answers is the following: + * Given a subdomain with associated DoFs, what is the largest subset of + * these DoFs that are allowed to be non-zero such that after calling + * ConstraintMatrix::distribute() the resulting solution vector will have + * support only within the given domain. Here @p cm is the ConstraintMatrix + * containing hanging nodes constraints. + * + * In case of parallel::distributed::Triangulation @p predicate will be called + * only for locally owned and ghost cells. The resulting index set may contain + * DoFs that are associated to the locally owned or ghost cells, but are not + * owned by the current MPI core. + */ + template + IndexSet + extract_dofs_with_support_on(const DoFHandlerType &dof_handler, + const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> &predicate, + const ConstraintMatrix &cm = ConstraintMatrix()); + /** * Extract a vector that represents the constant modes of the DoFHandler for * the components chosen by component_mask (see diff --git a/source/dofs/dof_tools.cc b/source/dofs/dof_tools.cc index 19314c16c1..796335ef7f 100644 --- a/source/dofs/dof_tools.cc +++ b/source/dofs/dof_tools.cc @@ -765,6 +765,110 @@ namespace DoFTools + template + IndexSet + extract_dofs_with_support_on(const DoFHandlerType &dof_handler, + const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> &predicate, + const ConstraintMatrix &cm) + { + const std::function< bool(const typename DoFHandlerType::active_cell_iterator &)> predicate_local + = [=] (const typename DoFHandlerType::active_cell_iterator & cell) -> bool + { + return cell->is_locally_owned() && predicate(cell); + }; + + std::vector local_dof_indices; + local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); + + // Get all the dofs that live on the subdomain: + std::set predicate_dofs; + + typename DoFHandlerType::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + if (!cell->is_artificial() && predicate(cell)) + { + local_dof_indices.resize (cell->get_fe().dofs_per_cell); + cell->get_dof_indices (local_dof_indices); + predicate_dofs.insert(local_dof_indices.begin(), + local_dof_indices.end()); + } + + // Get halo layer and accumulate its DoFs + std::set halo_dofs; + + const std::vector halo_cells = + GridTools::compute_active_cell_halo_layer(dof_handler, predicate_local); + for (typename std::vector::const_iterator it = halo_cells.begin(); + it != halo_cells.end(); ++it) + { + // skip halo cells that satisfy the given predicate and thereby + // shall be a part of the index set on another MPI core. + // Those could only be ghost cells with p::d::Tria. + if (predicate(*it)) + { + Assert ((*it)->is_ghost(), + ExcInternalError()); + continue; + } + + const unsigned int dofs_per_cell = (*it)->get_fe().dofs_per_cell; + local_dof_indices.resize (dofs_per_cell); + (*it)->get_dof_indices (local_dof_indices); + halo_dofs.insert(local_dof_indices.begin(), + local_dof_indices.end()); + } + + // A complication coming from the fact that DoFs living on locally + // owned cells which satisfy predicate may participate in constrains for + // DoFs outside of this region. + if (cm.n_constraints() > 0) + { + const std::vector > *line_ptr; + + // Remove DoFs that are part of constrains for DoFs outside + // of predicate. Since we will subtract halo_dofs from predicate_dofs, + // achieve this by extending halo_dofs with DoFs to which + // halo_dofs are constrained. + std::set extra_halo; + for (std::set::iterator it = halo_dofs.begin(); + it != halo_dofs.end(); ++it) + { + line_ptr = cm.get_constraint_entries(*it); + // if halo DoF is constrained, add all DoFs to which it's constrained. + if (line_ptr!=NULL) + { + const unsigned int line_size = line_ptr->size(); + for (unsigned int j=0; j &, const std::set &); + template + IndexSet + DoFTools::extract_dofs_with_support_on > + (const DoFHandler &, + const std::function< bool(const typename DoFHandler::active_cell_iterator &)> &, + const ConstraintMatrix&); + template void DoFTools::extract_hanging_node_dofs diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc new file mode 100644 index 0000000000..741020cd8c --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.cc @@ -0,0 +1,137 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +// + +// Test DoFTools::extract_dofs_with_support_on() + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +template +bool +pred_d(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49 && + cell->center()(1) < 0.49); +} + +template +bool +pred_r(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49 && + cell->center()(0) > 0.25 && + cell->center()(1) < 0.49); +} + + + +template +void test (const unsigned int flag) +{ + + // Setup system + Triangulation triangulation; + + GridGenerator::hyper_rectangle (triangulation, + Point(0,0), + Point(1,1)); + + if (flag == 0) + triangulation.refine_global(2); + else + triangulation.refine_global(1); + + // Extra refinement to generate hanging nodes + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if ( (flag==1 && pred_d(cell)) || + (flag==2 && !pred_d(cell)) ) + cell->set_refine_flag (); + + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement (); + + DoFHandler dh (triangulation); + + FE_Q fe(2); + dh.distribute_dofs (fe); + + ConstraintMatrix cm; + DoFTools::make_hanging_node_constraints(dh,cm); + + IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d, cm); + support.print(deallog); + + // print grid and DoFs for visual inspection + if (false) + { + std::cout << "-------------------- " << Utilities::int_to_string(dim) << Utilities::int_to_string(flag) << std::endl; + cm.print(std::cout); + + std::map > support_points; + MappingQ1 mapping; + DoFTools::map_dofs_to_support_points(mapping, dh, support_points); + + const std::string filename = + "grid" + Utilities::int_to_string(dim) + Utilities::int_to_string(flag) + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl + << "set output \"grid" << Utilities::int_to_string(dim) << Utilities::int_to_string(flag) << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" << std::endl; + GridOut().write_gnuplot (triangulation, f); + f << "e" << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, + support_points); + f << "e" << std::endl; + } + + dh.clear(); +} + +int main(int argc, char **argv) +{ + initlog(); + + test<2>(0); + test<2>(1); + test<2>(2); + + return 0; +} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output new file mode 100644 index 0000000000..d326d0283a --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_01.output @@ -0,0 +1,4 @@ + +DEAL::{[0,8], [12,14], [17,18], 20, 24} +DEAL::{[21,42]} +DEAL::{0, 4, 6, 8} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc new file mode 100644 index 0000000000..d74ba742fb --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.cc @@ -0,0 +1,185 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +// + +// similar to dof_tools_23, but for the p::d::Tria. As an MPI test, make sure +// that the total number of shape functions that are non-zero within the domain +// is the same. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +template +bool +pred_d(const typename parallel::distributed::Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.5 && + cell->center()(1) < 0.5); +} + + +std::string output_name(const unsigned int flag, const unsigned int subdomain) +{ + return "output" + Utilities::int_to_string(flag) + "_" + + Utilities::int_to_string(subdomain) + + ".vtu"; +} + + +template +void test (const unsigned int flag) +{ + // Setup system + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + + GridGenerator::hyper_rectangle (triangulation, + Point(0,0), + Point(1,1)); + + if (flag == 0) + triangulation.refine_global(2); + else + triangulation.refine_global(1); + + // Extra refinement to generate hanging nodes + for (typename parallel::distributed::Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned() && + ((flag==1 && pred_d(cell)) || + (flag==2 && !pred_d(cell))) ) + cell->set_refine_flag (); + + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement (); + + if (flag > 0) + triangulation.refine_global(1); + + DoFHandler dh (triangulation); + + FE_Q fe(2); + dh.distribute_dofs (fe); + + IndexSet locally_relevant_set; + DoFTools::extract_locally_relevant_dofs (dh, + locally_relevant_set); + + ConstraintMatrix cm; + cm.reinit (locally_relevant_set); + DoFTools::make_hanging_node_constraints(dh,cm); + cm.close (); + + const IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d, cm); + const IndexSet support_local = support & dh.locally_owned_dofs(); + + deallog << support.n_elements() << std::endl; + + // now accumulate the number of indices and make sure it's the same for + // various runs with different number of MPI cores + const unsigned int dofs_support = Utilities::MPI::sum(support_local.n_elements(), MPI_COMM_WORLD); + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + deallog << "accumulated: " << std::endl; + deallog << dofs_support << std::endl; + } + + // print grid and DoFs for visual inspection + if (false) + { + DataOut data_out; + data_out.attach_dof_handler (dh); + + std::vector > shape_functions(dh.n_dofs()); + for (unsigned int i = 0; i < dh.n_dofs(); ++i) + { + LinearAlgebra::distributed::Vector &s = shape_functions[i]; + s.reinit(dh.locally_owned_dofs(), + locally_relevant_set, + MPI_COMM_WORLD); + s = 0.; + if (dh.locally_owned_dofs().is_element(i)) + s[i] = 1.0; + s.compress(VectorOperation::insert); + cm.distribute(s); + s.update_ghost_values(); + + data_out.add_data_vector (s, + std::string("N_") + + Utilities::int_to_string(i)); + } + + Vector subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i filenames; + for (unsigned int i=0; i(0); + test<2>(1); + test<2>(2); + + deallog.pop(); +} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output new file mode 100644 index 0000000000..ad14a0a639 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=1.output @@ -0,0 +1,10 @@ + +DEAL:0:0::16 +DEAL:0:0::accumulated: +DEAL:0:0::16 +DEAL:0:0::76 +DEAL:0:0::accumulated: +DEAL:0:0::76 +DEAL:0:0::16 +DEAL:0:0::accumulated: +DEAL:0:0::16 diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output new file mode 100644 index 0000000000..24d5d3f17a --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=3.output @@ -0,0 +1,20 @@ + +DEAL:0:0::16 +DEAL:0:0::accumulated: +DEAL:0:0::16 +DEAL:0:0::61 +DEAL:0:0::accumulated: +DEAL:0:0::76 +DEAL:0:0::16 +DEAL:0:0::accumulated: +DEAL:0:0::16 + +DEAL:1:1::21 +DEAL:1:1::64 +DEAL:1:1::15 + + +DEAL:2:2::9 +DEAL:2:2::27 +DEAL:2:2::9 + diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output new file mode 100644 index 0000000000..4a18c5f509 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_02.mpirun=5.output @@ -0,0 +1,30 @@ + +DEAL:0:0::16 +DEAL:0:0::accumulated: +DEAL:0:0::16 +DEAL:0:0::49 +DEAL:0:0::accumulated: +DEAL:0:0::76 +DEAL:0:0::16 +DEAL:0:0::accumulated: +DEAL:0:0::16 + +DEAL:1:1::15 +DEAL:1:1::69 +DEAL:1:1::15 + + +DEAL:2:2::0 +DEAL:2:2::46 +DEAL:2:2::15 + + +DEAL:3:3::15 +DEAL:3:3::45 +DEAL:3:3::9 + + +DEAL:4:4::9 +DEAL:4:4::9 +DEAL:4:4::0 + diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc new file mode 100644 index 0000000000..e19e4ec036 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.cc @@ -0,0 +1,143 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +// + +// same as dof_tools_32 test DoFTools::extract_dofs_with_support_on() but +// for a slightly different configuration + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + + +template +bool +pred_left(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49); +} + +template +bool +pred_right(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) > 0.51); +} + + +template +bool +pred_r(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49 && + cell->center()(1) < 0.49) || + (cell->center()(0) > 0.49 && + cell->center()(1) > 0.49); +} + + + +template +void test (const bool left = true) +{ + + // Setup system + Triangulation triangulation; + + GridGenerator::hyper_rectangle (triangulation, + Point(0,0), + Point(1,1)); + + triangulation.refine_global(1); + + // Extra refinement to generate hanging nodes + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (pred_r(cell)) + cell->set_refine_flag (); + + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement (); + + DoFHandler dh (triangulation); + + FE_Q fe(2); + dh.distribute_dofs (fe); + + ConstraintMatrix cm; + DoFTools::make_hanging_node_constraints(dh,cm); + + IndexSet support = left ? + DoFTools::extract_dofs_with_support_on(dh, pred_left, cm) : + DoFTools::extract_dofs_with_support_on(dh, pred_right, cm); + support.print(deallog); + + // print grid and DoFs for visual inspection + if (false) + { + std::cout << "-------------------- " << Utilities::int_to_string(dim) << std::endl; + cm.print(std::cout); + + std::map > support_points; + MappingQ1 mapping; + DoFTools::map_dofs_to_support_points(mapping, dh, support_points); + + const std::string filename = + "grid" + Utilities::int_to_string(dim) + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl + << "set output \"grid" << Utilities::int_to_string(dim) << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" << std::endl; + GridOut().write_gnuplot (triangulation, f); + f << "e" << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, + support_points); + f << "e" << std::endl; + } + + dh.clear(); +} + +int main(int argc, char **argv) +{ + initlog(); + + test<2>(true); + test<2>(false); + + return 0; +} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output new file mode 100644 index 0000000000..76ee14d0e0 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_03.output @@ -0,0 +1,3 @@ + +DEAL::{[9,10], 12, [14,38]} +DEAL::{1, 3, [5,8], [39,60]} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc new file mode 100644 index 0000000000..7126ab0558 --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.cc @@ -0,0 +1,327 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- +// + +// test DoFTools::extract_dofs_with_support_on() for calculation of the RHS +// function in parallel + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + + +std::string output_name(const unsigned int subdomain) +{ + return "output_" + + Utilities::int_to_string(subdomain) + + ".vtu"; +} + +template +bool +pred_d(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49); +} + +template +bool +pred_r(const typename Triangulation::active_cell_iterator &cell) +{ + return (cell->center()(0) < 0.49 && + cell->center()(1) < 0.49) || + (cell->center()(0) > 0.49 && + cell->center()(1) > 0.49); +} + + + +template +void test () +{ + + // Setup system + parallel::distributed::Triangulation triangulation(MPI_COMM_WORLD); + + GridGenerator::hyper_rectangle (triangulation, + Point(0,0), + Point(1,1)); + + triangulation.refine_global(1); + + // Extra refinement to generate hanging nodes + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + if (cell->is_locally_owned() && pred_r(cell)) + cell->set_refine_flag (); + + triangulation.prepare_coarsening_and_refinement(); + triangulation.execute_coarsening_and_refinement (); + + DoFHandler dh (triangulation); + + FE_Q fe(2); + dh.distribute_dofs (fe); + + const IndexSet &locally_owned_set = dh.locally_owned_dofs(); + IndexSet locally_relevant_set; + DoFTools::extract_locally_relevant_dofs (dh, + locally_relevant_set); + + ConstraintMatrix cm; + cm.reinit (locally_relevant_set); + DoFTools::make_hanging_node_constraints(dh,cm); + cm.close (); + + std::vector local_dof_indices(fe.dofs_per_cell); + + // get support on the predicate + IndexSet support = DoFTools::extract_dofs_with_support_on(dh, pred_d, cm); + IndexSet local_support = support & locally_owned_set; + + // rhs vectors: + LinearAlgebra::distributed::Vector sparse_rhs; + LinearAlgebra::distributed::Vector rhs; + sparse_rhs.reinit(locally_owned_set, locally_relevant_set, MPI_COMM_WORLD); + rhs.reinit( locally_owned_set, locally_relevant_set, MPI_COMM_WORLD); + + sparse_rhs.zero_out_ghosts(); + rhs.zero_out_ghosts(); + + // assemble RHS which has a local support: + const std::function< double(const Point &)> rhs_func + = [=] (const Point &p) -> double + { + return p[0] > 0.5 ? 0. : 0.5 - p[0]; + }; + + Vector local_rhs(fe.dofs_per_cell); + QGauss quadrature(3); + FEValues fe_values(fe, quadrature, + update_values | update_JxW_values | update_quadrature_points); + for (typename DoFHandler::active_cell_iterator + cell = dh.begin_active(); + cell != dh.end(); ++cell) + if (cell->is_locally_owned() && pred_d(cell)) + { + fe_values.reinit(cell); + cell->get_dof_indices (local_dof_indices); + + const std::vector> &q_points = fe_values.get_quadrature_points(); + + local_rhs = 0.; + for (unsigned int i = 0; i < fe.dofs_per_cell; ++i) + for (unsigned int q = 0; q < quadrature.size(); ++q) + local_rhs[i] += fe_values.shape_value (i, q) * + rhs_func(q_points[q]) * + fe_values.JxW (q); + + cm.distribute_local_to_global (local_rhs, + local_dof_indices, + rhs); + + // copy-paste of CM distribute_local_to_global and + // add is_element() checks: + auto local_vector_begin = local_rhs.begin(); + const auto local_vector_end = local_rhs.end(); + auto local_indices_begin = local_dof_indices.begin(); + const std::vector > *line_ptr; + for ( ; local_vector_begin != local_vector_end; + ++local_vector_begin, ++local_indices_begin) + { + line_ptr = cm.get_constraint_entries(*local_indices_begin); + if (line_ptr==NULL) // unconstrained + { + if (support.is_element(*local_indices_begin)) + sparse_rhs(*local_indices_begin) += *local_vector_begin; + } + else + { + const unsigned int line_size = line_ptr->size(); + for (unsigned int j=0; j > support_points; + MappingQ1 mapping; + DoFTools::map_dofs_to_support_points(mapping, dh, support_points); + + const std::string filename = + "grid" + dealii::Utilities::int_to_string(n_mpi_processes) + dealii::Utilities::int_to_string(this_mpi_process); + std::ofstream f((filename+".gp").c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" << std::endl + << "set output \"" << filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" << std::endl; + GridOut().write_gnuplot (triangulation, f); + f << "e" << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, + support_points); + f << "e" << std::endl; + + } + + { + DataOut data_out; + data_out.attach_dof_handler (dh); + + std::vector > shape_functions(dh.n_dofs()); + for (unsigned int i = 0; i < dh.n_dofs(); ++i) + { + LinearAlgebra::distributed::Vector sl(locally_owned_set, MPI_COMM_WORLD); + sl = 0.; + if (locally_owned_set.is_element(i)) + sl[i] = 1.0; + cm.distribute(sl); + + LinearAlgebra::distributed::Vector &s = shape_functions[i]; + s.reinit(locally_owned_set, + locally_relevant_set, + MPI_COMM_WORLD); + s = 0.; + s = sl; + + data_out.add_data_vector (s, + std::string("N_") + + dealii::Utilities::int_to_string(i)); + } + + Vector subdomain (triangulation.n_active_cells()); + for (unsigned int i=0; i filenames; + for (unsigned int i=0; i(); + + return 0; +} diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output new file mode 100644 index 0000000000..04c93dbf6b --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.mpirun=2.output @@ -0,0 +1 @@ +DEAL::Ok diff --git a/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output new file mode 100644 index 0000000000..04c93dbf6b --- /dev/null +++ b/tests/dofs/dof_tools_extract_dofs_with_support_on_04.output @@ -0,0 +1 @@ +DEAL::Ok -- 2.39.5