Line data Source code
1 : /*
2 : * Copyright (c) 2012-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__FUNCTION_SPACES__INTEGRATE_FLUX__
34 : #define __H__UG__LIB_DISC__FUNCTION_SPACES__INTEGRATE_FLUX__
35 :
36 : #include <cmath>
37 :
38 : #include "common/common.h"
39 : #include "lib_grid/tools/subset_group.h"
40 : #include "lib_disc/common/function_group.h"
41 : #include "lib_disc/common/groups_util.h"
42 :
43 : #include "lib_disc/assemble_interface.h"
44 : #include "lib_disc/spatial_disc/elem_disc/elem_disc_interface.h"
45 : #include "lib_disc/spatial_disc/constraints/constraint_interface.h"
46 :
47 : namespace ug{
48 :
49 : template <typename TGridFunction, int dim>
50 0 : void IntegrateDiscFlux(std::vector<number>& vValue,
51 : const TGridFunction& rDefect,
52 : const FunctionGroup& fctGrp,
53 : int si)
54 : {
55 : // get iterators for all elems on subset
56 : typedef typename TGridFunction::template dim_traits<dim>::const_iterator const_iterator;
57 : typedef typename TGridFunction::template dim_traits<dim>::grid_base_object element_type;
58 : const_iterator iter = rDefect.template begin<element_type>(si);
59 : const_iterator iterEnd = rDefect.template end<element_type>(si);
60 :
61 : // return if no dofs on this element types
62 0 : if(rDefect.max_dofs(dim) == 0) return;
63 :
64 : // check sizes
65 0 : if(vValue.size() != fctGrp.size())
66 0 : UG_THROW("IntegrateDiscFlux: number of values and number of function mismatch.");
67 :
68 : // allocate memory for indices
69 0 : std::vector<std::vector<DoFIndex> > vInd(fctGrp.size());
70 :
71 : // loop elements of subset
72 0 : for( ; iter != iterEnd; ++iter)
73 : {
74 : // get element
75 : element_type* elem = *iter;
76 :
77 : // get multi-indices of elem
78 0 : for(size_t f = 0; f < fctGrp.size(); ++f)
79 : rDefect.inner_dof_indices(elem, fctGrp[f], vInd[f]);
80 :
81 : // sum values
82 0 : for(size_t f = 0; f < fctGrp.size(); ++f)
83 0 : for(size_t i = 0; i < vInd[f].size(); ++i)
84 0 : vValue[f] += DoFRef( rDefect, vInd[f][i]);
85 : }
86 0 : }
87 :
88 :
89 : template <typename TGridFunction>
90 0 : number IntegrateDiscFlux(SmartPtr<IAssemble<typename TGridFunction::algebra_type> > spAssemble,
91 : SmartPtr<TGridFunction> spGridFct,
92 : const char* pCmps,
93 : const char* subsets)
94 : {
95 : // read subsets
96 0 : SubsetGroup ssGrp(spGridFct->domain()->subset_handler());
97 0 : if(subsets != NULL)
98 : {
99 0 : ssGrp.add(TokenizeString(subsets));
100 0 : if(!SameDimensionsInAllSubsets(ssGrp))
101 0 : UG_THROW("IntegrateDiscFlux: Subsets '"<<subsets<<"' do not have same dimension."
102 : "Can not integrate on subsets of different dimensions.");
103 : }
104 : else
105 : {
106 : // add all subsets and remove lower dim subsets afterwards
107 0 : ssGrp.add_all();
108 0 : RemoveLowerDimSubsets(ssGrp);
109 : }
110 :
111 : // read functions
112 0 : FunctionGroup fctGrp(spGridFct->function_pattern(), TokenizeString(pCmps));
113 :
114 : // check if something to do
115 0 : if(fctGrp.size() == 0) return 0;
116 :
117 : // create vector
118 0 : TGridFunction Defect(*spGridFct);
119 :
120 : // remember enabled-flags
121 : // get assemble adapter
122 0 : SmartPtr<AssemblingTuner<typename TGridFunction::algebra_type> > spAssTuner = spAssemble->ass_tuner();
123 : const int ElemTypesEnabled = spAssTuner->enabled_elem_discs();
124 : const int ConstraintTypesEnabled = spAssTuner->enabled_constraints();
125 :
126 : // remove bnd components
127 0 : spAssTuner->enable_elem_discs(ElemTypesEnabled & (~EDT_BND));
128 0 : spAssTuner->enable_constraints(ConstraintTypesEnabled & (~CT_DIRICHLET));
129 :
130 : // compute defect
131 0 : spAssemble->assemble_defect(Defect, *spGridFct);
132 :
133 : // reset flags to original status
134 : spAssTuner->enable_elem_discs(ElemTypesEnabled);
135 : spAssTuner->enable_constraints(ConstraintTypesEnabled);
136 :
137 : // reset value
138 0 : std::vector<number> vValue(fctGrp.size(), 0);
139 :
140 : // loop subsets
141 0 : for(size_t i = 0; i < ssGrp.size(); ++i)
142 : {
143 : // get subset index
144 : const int si = ssGrp[i];
145 :
146 : // error if function is not defined in subset
147 0 : for(size_t f = 0; f < fctGrp.size(); ++f)
148 : {
149 0 : if(!spGridFct->is_def_in_subset(fctGrp[f], si))
150 0 : UG_THROW("IntegrateDiscFlux: Function "<<fctGrp[f]<<" not defined "
151 : "on subset "<<si<<".");
152 : }
153 :
154 : // check dimension
155 0 : if(ssGrp.dim(i) > TGridFunction::dim)
156 0 : UG_THROW("IntegrateDiscFlux: Dimension of subset is "<<ssGrp.dim(i)<<", but "
157 : " World Dimension is "<<TGridFunction::dim<<". Cannot integrate this.");
158 :
159 : // integrate elements of subset
160 : try{
161 0 : for(int d = 0; d < ssGrp.dim(i)+1; ++d)
162 : {
163 0 : if(d==0) IntegrateDiscFlux<TGridFunction, 0>(vValue, Defect, fctGrp, si);
164 0 : else if(d==1) IntegrateDiscFlux<TGridFunction, 1>(vValue, Defect, fctGrp, si);
165 0 : else if(d==2) IntegrateDiscFlux<TGridFunction, 2>(vValue, Defect, fctGrp, si);
166 0 : else if(d==3) IntegrateDiscFlux<TGridFunction, 3>(vValue, Defect, fctGrp, si);
167 0 : else UG_THROW("IntegrateDiscFlux: Dimension "<<d<<" not supported. "
168 : " World dimension is "<<TGridFunction::dim<<".");
169 : }
170 : }
171 0 : UG_CATCH_THROW("IntegrateSubsets: Integration failed on subset "<<si);
172 :
173 : }
174 :
175 : #ifdef UG_PARALLEL
176 : // sum over processes
177 : if(pcl::NumProcs() > 1)
178 : {
179 : pcl::ProcessCommunicator com;
180 : for(size_t f = 0; f < vValue.size(); ++f)
181 : {
182 : number local = vValue[f];
183 : com.allreduce(&local, &vValue[f], 1, PCL_DT_DOUBLE, PCL_RO_SUM);
184 : }
185 : }
186 : #endif
187 :
188 : // return the result
189 0 : if(vValue.size() == 1) return vValue[0];
190 : else
191 : {
192 : number sum = 0;
193 0 : for(size_t f = 0; f < vValue.size(); ++f){
194 0 : UG_LOG("Integral for '"<<fctGrp.name(f)<<"': "<<vValue[f]<<"\n");
195 0 : sum += vValue[f];
196 : }
197 0 : return sum;
198 : }
199 0 : }
200 :
201 : } // end namespace ug
202 :
203 : #endif /* __H__UG__LIB_DISC__FUNCTION_SPACES__INTEGRATE_FLUX__ */
|
