Line data Source code
1 : /*
2 : * Copyright (c) 2013-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 : #include "common/util/string_util.h"
34 : #include "quadrature_provider.h"
35 : #include "gauss/gauss_quad.h"
36 : #include "gauss/gauss_quad_vertex.h"
37 : #include "newton_cotes/newton_cotes.h"
38 : #include "gauss_legendre/gauss_legendre.h"
39 : #include "gauss_jacobi/gauss_jacobi10.h"
40 : #include "gauss_jacobi/gauss_jacobi20.h"
41 : #include "gauss_tensor_prod/gauss_tensor_prod.h"
42 : #include "lib_disc/reference_element/reference_element.h"
43 : #include <algorithm>
44 : #include <locale>
45 :
46 : namespace ug{
47 :
48 : ////////////////////////////////////////////////////////////////////////////////
49 : // gauss
50 : ////////////////////////////////////////////////////////////////////////////////
51 :
52 : template <>
53 : const QuadratureRule<0>*
54 0 : QuadratureRuleProvider<0>::create_gauss_rule(ReferenceObjectID roid,
55 : size_t order)
56 : {
57 : QuadratureRule<0>* q = NULL;
58 : try{
59 0 : switch(roid){
60 0 : case ROID_VERTEX: q = new GaussQuadratureVertex(); break;
61 0 : default: UG_THROW("QuadratureRuleProvider<"<<dim<<">: "<<roid<<" not supported.");
62 : }
63 0 : }catch(...){return NULL;}
64 0 : return q;
65 : }
66 :
67 : template <>
68 : const QuadratureRule<1>*
69 0 : QuadratureRuleProvider<1>::create_gauss_rule(ReferenceObjectID roid,
70 : size_t order)
71 : {
72 : QuadratureRule<1>* q = NULL;
73 : try{
74 0 : switch(roid){
75 0 : case ROID_EDGE: q = new FlexGaussQuadrature<ReferenceEdge>(order); break;
76 0 : default: UG_THROW("QuadratureRuleProvider<"<<dim<<">: "<<roid<<" not supported.");
77 : }
78 0 : }catch(...){return NULL;}
79 : return q;
80 : }
81 :
82 : template <>
83 : const QuadratureRule<2>*
84 0 : QuadratureRuleProvider<2>::create_gauss_rule(ReferenceObjectID roid,
85 : size_t order)
86 : {
87 : QuadratureRule<2>* q = NULL;
88 : try{
89 0 : switch(roid){
90 0 : case ROID_TRIANGLE: q = new FlexGaussQuadrature<ReferenceTriangle>(order); break;
91 0 : case ROID_QUADRILATERAL: q = new FlexGaussQuadrature<ReferenceQuadrilateral>(order); break;
92 0 : default: UG_THROW("QuadratureRuleProvider<"<<dim<<">: "<<roid<<" not supported.");
93 : }
94 0 : }catch(...){return NULL;}
95 : return q;
96 : }
97 :
98 : template <>
99 : const QuadratureRule<3>*
100 0 : QuadratureRuleProvider<3>::create_gauss_rule(ReferenceObjectID roid,
101 : size_t order)
102 : {
103 : QuadratureRule<3>* q = NULL;
104 : try{
105 0 : switch(roid){
106 0 : case ROID_TETRAHEDRON: q = new FlexGaussQuadrature<ReferenceTetrahedron>(order); break;
107 0 : case ROID_PYRAMID: q = new FlexGaussQuadrature<ReferencePyramid>(order); break;
108 0 : case ROID_PRISM: q = new FlexGaussQuadrature<ReferencePrism>(order); break;
109 0 : case ROID_HEXAHEDRON: q = new FlexGaussQuadrature<ReferenceHexahedron>(order); break;
110 0 : case ROID_OCTAHEDRON: q = new FlexGaussQuadrature<ReferenceOctahedron>(order); break;
111 0 : default: UG_THROW("QuadratureRuleProvider<"<<dim<<">: "<<roid<<" not supported.");
112 : }
113 0 : }catch(...){return NULL;}
114 : return q;
115 : }
116 :
117 : ////////////////////////////////////////////////////////////////////////////////
118 : // gauss-legendre
119 : ////////////////////////////////////////////////////////////////////////////////
120 :
121 : template <>
122 : const QuadratureRule<0>*
123 0 : QuadratureRuleProvider<0>::create_gauss_legendre_rule(ReferenceObjectID roid, size_t order)
124 : {
125 0 : return NULL;
126 : }
127 :
128 : template <>
129 : const QuadratureRule<1>*
130 0 : QuadratureRuleProvider<1>::create_gauss_legendre_rule(ReferenceObjectID roid, size_t order)
131 : {
132 : QuadratureRule<1>* q = NULL;
133 : try{
134 0 : q = new GaussLegendre(order);
135 0 : }catch(...){return NULL;}
136 : return q;
137 : }
138 :
139 : template <>
140 : const QuadratureRule<2>*
141 0 : QuadratureRuleProvider<2>::create_gauss_legendre_rule(ReferenceObjectID roid,
142 : size_t order)
143 : {
144 : QuadratureRule<2>* q = NULL;
145 : try{
146 0 : switch(roid){
147 0 : case ROID_QUADRILATERAL: q = new GaussQuadratureQuadrilateral(order); break;
148 0 : case ROID_TRIANGLE: q = new GaussQuadratureTriangle(order); break;
149 0 : default: UG_THROW("QuadratureRuleProvider<"<<dim<<">: "<<roid<<" not supported.");
150 : }
151 0 : }catch(...){return NULL;}
152 : return q;
153 : }
154 :
155 : template <>
156 : const QuadratureRule<3>*
157 0 : QuadratureRuleProvider<3>::create_gauss_legendre_rule(ReferenceObjectID roid,
158 : size_t order)
159 : {
160 : QuadratureRule<3>* q = NULL;
161 : try{
162 0 : switch(roid){
163 0 : case ROID_TETRAHEDRON: q = new GaussQuadratureTetrahedron(order); break;
164 0 : case ROID_PRISM: q = new GaussQuadraturePrism(order); break;
165 0 : case ROID_PYRAMID: q = new GaussQuadraturePyramid(order); break;
166 0 : case ROID_HEXAHEDRON: q = new GaussQuadratureHexahedron(order); break;
167 0 : case ROID_OCTAHEDRON: q = new GaussQuadratureOctahedron(order); break;
168 0 : default: UG_THROW("QuadratureRuleProvider<"<<dim<<">: "<<roid<<" not supported.");
169 : }
170 0 : }catch(...){return NULL;}
171 : return q;
172 : }
173 :
174 : ////////////////////////////////////////////////////////////////////////////////
175 : // newton-cotes
176 : ////////////////////////////////////////////////////////////////////////////////
177 :
178 : template <>
179 : const QuadratureRule<1>*
180 0 : QuadratureRuleProvider<1>::create_newton_cotes_rule(ReferenceObjectID roid, size_t order)
181 : {
182 : QuadratureRule<1>* q = NULL;
183 : try{
184 0 : q = new NewtonCotes(order);
185 0 : }catch(...){return NULL;}
186 : return q;
187 : }
188 :
189 : template <int TDim>
190 : const QuadratureRule<TDim>*
191 0 : QuadratureRuleProvider<TDim>::create_newton_cotes_rule(ReferenceObjectID roid, size_t order)
192 : {
193 0 : return NULL;
194 : }
195 :
196 : ////////////////////////////////////////////////////////////////////////////////
197 : // general
198 : ////////////////////////////////////////////////////////////////////////////////
199 :
200 : template <int TDim>
201 0 : QuadratureRuleProvider<TDim>::QuadratureRuleProvider()
202 : {
203 0 : for(int type = 0; type < NUM_QUADRATURE_TYPES; ++type)
204 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid)
205 : m_vRule[type][roid].clear();
206 0 : }
207 :
208 : template <int TDim>
209 0 : QuadratureRuleProvider<TDim>::~QuadratureRuleProvider()
210 : {
211 0 : for(int type = 0; type < NUM_QUADRATURE_TYPES; ++type)
212 0 : for(int roid = 0; roid < NUM_REFERENCE_OBJECTS; ++roid)
213 0 : for(size_t order = 0; order < m_vRule[type][roid].size(); ++order)
214 0 : if(m_vRule[type][roid][order] != NULL)
215 0 : delete m_vRule[type][roid][order];
216 0 : }
217 :
218 : template <int TDim>
219 : const QuadratureRule<TDim>&
220 0 : QuadratureRuleProvider<TDim>::get_quad_rule(ReferenceObjectID roid,
221 : size_t order,
222 : QuadType type)
223 : {
224 : // check if order present, else resize and create
225 0 : if(order >= m_vRule[type][roid].size() ||
226 0 : m_vRule[type][roid][order] == NULL)
227 0 : create_rule(roid, order, type);
228 :
229 : // return correct order
230 0 : return *m_vRule[type][roid][order];
231 : }
232 :
233 : template <int TDim>
234 : void
235 0 : QuadratureRuleProvider<TDim>::create_rule(ReferenceObjectID roid,
236 : size_t order,
237 : QuadType type)
238 : {
239 : // resize vector if needed
240 0 : if(m_vRule[type][roid].size() <= order) m_vRule[type][roid].resize(order+1, NULL);
241 0 : if(m_vRule[type][roid][order] != NULL)
242 0 : delete m_vRule[type][roid][order];
243 0 : m_vRule[type][roid][order] = NULL;
244 :
245 0 : switch(type){
246 0 : case BEST: {
247 : // 1. Try GaussQuad
248 0 : m_vRule[type][roid][order] = create_gauss_rule(roid, order);
249 0 : if(m_vRule[type][roid][order] != NULL) break;
250 :
251 : // 2. Try Newton-Cotes
252 0 : m_vRule[type][roid][order] = create_newton_cotes_rule(roid, order);
253 0 : if(m_vRule[type][roid][order] != NULL) break;
254 :
255 : // 3. Try Gauss-Legendre
256 0 : m_vRule[type][roid][order] = create_gauss_legendre_rule(roid, order);
257 : if(m_vRule[type][roid][order] != NULL) break;
258 : }break;
259 0 : case GAUSS: {
260 0 : m_vRule[type][roid][order] = create_gauss_rule(roid, order);
261 0 : }break;
262 0 : case GAUSS_LEGENDRE: {
263 0 : m_vRule[type][roid][order] = create_gauss_legendre_rule(roid, order);
264 0 : }break;
265 0 : case NEWTON_COTES: {
266 0 : m_vRule[type][roid][order] = create_newton_cotes_rule(roid, order);
267 0 : }break;
268 0 : default: UG_THROW("QuadratureRuleProvider<"<<dim<<">: Cannot create rule for "
269 : <<roid<<", order "<<order<<" and type "<<type);
270 : }
271 :
272 0 : if(m_vRule[type][roid][order] == NULL)
273 0 : UG_THROW("QuadratureRuleProvider<"<<dim<<">: Cannot create rule for "
274 : <<roid<<", order "<<order<<" and type "<<type);
275 0 : }
276 :
277 : template <int TDim>
278 : const QuadratureRule<TDim>&
279 0 : QuadratureRuleProvider<TDim>::get(ReferenceObjectID roid, size_t order,
280 : QuadType type)
281 : {
282 : // forward request
283 0 : return instance().get_quad_rule(roid, order, type);
284 : }
285 :
286 0 : std::ostream& operator<<(std::ostream& out, const QuadType& v)
287 : {
288 0 : switch(v)
289 : {
290 0 : case BEST: out << "Best"; break;
291 0 : case GAUSS: out << "Gauss"; break;
292 0 : case NEWTON_COTES: out << "Newton-Cotes"; break;
293 0 : case GAUSS_LEGENDRE: out << "Gauss-Legendre"; break;
294 0 : default: out << "invalid";
295 : }
296 0 : return out;
297 : }
298 :
299 0 : QuadType GetQuadratureType(const std::string& name)
300 : {
301 0 : std::string n = TrimString(name);
302 : std::transform(n.begin(), n.end(), n.begin(), ::tolower);
303 0 : if(n == "best") return BEST;
304 0 : if(n == "gauss") return GAUSS;
305 0 : if(n == "gauss-legendre") return GAUSS_LEGENDRE;
306 0 : if(n == "newton-cotes") return NEWTON_COTES;
307 :
308 0 : UG_THROW("GetQuadratureType: type '"<<name<<"' not recognized. Options "
309 : "are: best, gauss, gauss-legendre, newton-cotes.");
310 : }
311 :
312 :
313 : ////////////////////////////////////////////////////////////////////////////////
314 : // explicit template instantiations
315 : ////////////////////////////////////////////////////////////////////////////////
316 :
317 : template class QuadratureRuleProvider<0>;
318 : template class QuadratureRuleProvider<1>;
319 : template class QuadratureRuleProvider<2>;
320 : template class QuadratureRuleProvider<3>;
321 :
322 : } // namespace ug
|
