Line data Source code
1 : /*
2 : * Copyright (c) 2013-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Lisa Grau
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 :
34 : #include "../quadrature.h"
35 : #include "common/util/provider.h"
36 : #include "gauss_tensor_prod.h"
37 : #include "../gauss_legendre/gauss_legendre.h"
38 : #include "../gauss_jacobi/gauss_jacobi10.h"
39 : #include "../gauss_jacobi/gauss_jacobi20.h"
40 : #include "lib_disc/reference_element/reference_mapping_provider.h"
41 :
42 : namespace ug {
43 :
44 0 : GaussQuadratureTriangle::GaussQuadratureTriangle(size_t order)
45 : {
46 0 : GaussLegendre quadRule = GaussLegendre(order);
47 0 : GaussJacobi10 quadRule10 = GaussJacobi10(order);
48 :
49 0 : m_order = std::min(quadRule.order(), quadRule10.order());
50 0 : m_numPoints = quadRule10.size() * quadRule.size();
51 0 : position_type* pvPoint = new position_type[m_numPoints];
52 0 : weight_type* pvWeight = new weight_type[m_numPoints];
53 :
54 : size_t cnt = 0;
55 0 : for(size_t i = 0; i < quadRule.size(); i++){
56 0 : for(size_t j = 0; j < quadRule10.size(); j++, cnt++){
57 0 : pvPoint[cnt][0] = quadRule10.point(j)[0];
58 0 : pvPoint[cnt][1] = (1.0 - quadRule10.point(j)[0] ) * quadRule.point(i)[0];
59 0 : pvWeight[cnt] = quadRule.weight(i) * quadRule10.weight(j);
60 : }
61 : }
62 0 : m_pvPoint = pvPoint;
63 0 : m_pvWeight = pvWeight;
64 0 : };
65 :
66 0 : GaussQuadratureTriangle::~GaussQuadratureTriangle()
67 : {
68 0 : delete[] m_pvPoint;
69 0 : delete[] m_pvWeight;
70 0 : }
71 :
72 0 : GaussQuadratureQuadrilateral::GaussQuadratureQuadrilateral(size_t order)
73 : {
74 0 : GaussLegendre quadRule = GaussLegendre(order);
75 :
76 0 : m_order = std::min(quadRule.order(), quadRule.order());
77 0 : m_numPoints = quadRule.size() * quadRule.size();
78 0 : position_type* pvPoint = new position_type[m_numPoints];
79 0 : weight_type* pvWeight = new weight_type[m_numPoints];
80 :
81 : size_t cnt = 0;
82 0 : for(size_t i = 0; i < quadRule.size(); i ++) {
83 0 : for(size_t j = 0; j < quadRule.size(); j++, cnt++) {
84 0 : pvPoint[cnt][0] = quadRule.point(i)[0];
85 0 : pvPoint[cnt][1] = quadRule.point(j)[0];
86 0 : pvWeight[cnt] = quadRule.weight(i) * quadRule.weight(j);
87 : }
88 : }
89 0 : m_pvPoint = pvPoint;
90 0 : m_pvWeight = pvWeight;
91 0 : };
92 :
93 0 : GaussQuadratureQuadrilateral::~GaussQuadratureQuadrilateral()
94 : {
95 0 : delete[] m_pvPoint;
96 0 : delete[] m_pvWeight;
97 0 : }
98 :
99 0 : GaussQuadratureHexahedron::GaussQuadratureHexahedron(size_t order)
100 : {
101 0 : GaussLegendre quadRule = GaussLegendre(order);
102 :
103 0 : m_order = std::min(quadRule.order(), std::min(quadRule.order(), quadRule.order()));
104 0 : m_numPoints = quadRule.size() * quadRule.size() * quadRule.size();
105 0 : position_type* pvPoint = new position_type[m_numPoints];
106 0 : weight_type* pvWeight = new weight_type[m_numPoints];
107 :
108 : size_t cnt = 0;
109 0 : for(size_t i = 0; i < quadRule.size(); i ++) {
110 0 : for(size_t j = 0; j < quadRule.size(); j++) {
111 0 : for(size_t k = 0; k < quadRule.size(); k++, cnt++) {
112 0 : pvPoint[cnt][0] = quadRule.point(i)[0];
113 0 : pvPoint[cnt][1] = quadRule.point(j)[0];
114 0 : pvPoint[cnt][2] = quadRule.point(k)[0];
115 0 : pvWeight[cnt] = quadRule.weight(i) * quadRule.weight(j) * quadRule.weight(k);
116 : }
117 : }
118 : }
119 0 : m_pvPoint = pvPoint;
120 0 : m_pvWeight = pvWeight;
121 0 : };
122 :
123 0 : GaussQuadratureHexahedron::~GaussQuadratureHexahedron()
124 : {
125 0 : delete[] m_pvPoint;
126 0 : delete[] m_pvWeight;
127 0 : }
128 :
129 0 : GaussQuadratureTetrahedron::GaussQuadratureTetrahedron(size_t order)
130 : {
131 0 : GaussLegendre quadRule = GaussLegendre(order);
132 0 : GaussJacobi10 quadRule10 = GaussJacobi10(order);
133 0 : GaussJacobi20 quadRule20 = GaussJacobi20(order);
134 :
135 0 : m_order = std::min(quadRule.order(), std::min(quadRule10.order(), quadRule20.order()));
136 0 : m_numPoints = quadRule.size() * quadRule10.size() * quadRule20.size();
137 0 : position_type* pvPoint = new position_type[m_numPoints];
138 0 : weight_type* pvWeight = new weight_type[m_numPoints];
139 :
140 : size_t cnt = 0;
141 0 : for(size_t i = 0; i < quadRule20.size(); i++) {
142 0 : for(size_t j = 0; j < quadRule10.size(); j++) {
143 0 : for(size_t k = 0; k < quadRule.size(); k++, cnt++) {
144 0 : pvPoint[cnt][0] = quadRule20.point(i)[0];
145 0 : pvPoint[cnt][1] = (1.0 - quadRule20.point(i)[0] ) * quadRule10.point(j)[0];
146 0 : pvPoint[cnt][2] = (1.0 - quadRule20.point(i)[0]) * (1.0 - quadRule10.point(j)[0]) * quadRule.point(k)[0];
147 0 : pvWeight[cnt] = quadRule20.weight(i) * quadRule10.weight(j) * quadRule.weight(k);
148 : }
149 : }
150 : }
151 0 : m_pvPoint = pvPoint;
152 0 : m_pvWeight = pvWeight;
153 0 : };
154 :
155 0 : GaussQuadratureTetrahedron::~GaussQuadratureTetrahedron()
156 : {
157 0 : delete[] m_pvPoint;
158 0 : delete[] m_pvWeight;
159 0 : }
160 :
161 0 : GaussQuadraturePrism::GaussQuadraturePrism(size_t order)
162 : {
163 0 : GaussLegendre quadRule = GaussLegendre(order);
164 0 : GaussJacobi10 quadRule10 = GaussJacobi10(order);
165 :
166 0 : m_order = std::min(quadRule.order(), quadRule10.order());
167 0 : m_numPoints = quadRule10.size() * quadRule.size() * quadRule.size();
168 0 : position_type* pvPoint = new position_type[m_numPoints];
169 0 : weight_type* pvWeight = new weight_type[m_numPoints];
170 :
171 : size_t cnt = 0;
172 0 : for(size_t i = 0; i < quadRule10.size(); i++) {
173 0 : for(size_t j = 0; j < quadRule.size(); j++) {
174 0 : for(size_t k = 0; k < quadRule.size(); k++, cnt++) {
175 0 : pvPoint[cnt][0] = quadRule10.point(i)[0];
176 0 : pvPoint[cnt][1] = (1.0 - quadRule10.point(i)[0]) * quadRule.point(j)[0];
177 0 : pvPoint[cnt][2] = quadRule.point(k)[0];
178 0 : pvWeight[cnt] = quadRule10.weight(i) * quadRule.weight(j) * quadRule.weight(k);
179 : }
180 : }
181 : }
182 0 : m_pvPoint = pvPoint;
183 0 : m_pvWeight = pvWeight;
184 0 : };
185 :
186 0 : GaussQuadraturePrism::~GaussQuadraturePrism()
187 : {
188 0 : delete[] m_pvPoint;
189 0 : delete[] m_pvWeight;
190 0 : }
191 :
192 0 : GaussQuadraturePyramid::GaussQuadraturePyramid(size_t order)
193 : {
194 0 : GaussQuadratureTetrahedron quadRule = GaussQuadratureTetrahedron(order);
195 :
196 0 : m_order = quadRule.order();
197 0 : m_numPoints = quadRule.size() * 2;
198 0 : position_type* pvPoint = new position_type[m_numPoints];
199 0 : weight_type* pvWeight = new weight_type[m_numPoints];
200 :
201 0 : MathVector<3> Tet1Co[4];
202 : Tet1Co[0] = MathVector<3>(0,0,0);
203 : Tet1Co[1] = MathVector<3>(1,1,0);
204 : Tet1Co[2] = MathVector<3>(0,0,1);
205 : Tet1Co[3] = MathVector<3>(0,1,0);
206 :
207 0 : MathVector<3> Tet2Co[4];
208 : Tet2Co[0] = MathVector<3>(0,0,0);
209 : Tet2Co[1] = MathVector<3>(1,0,0);
210 : Tet2Co[2] = MathVector<3>(1,1,0);
211 : Tet2Co[3] = MathVector<3>(0,0,1);
212 :
213 : DimReferenceMapping<3, 3>& map1 =
214 : ReferenceMappingProvider::get<3,3>(ROID_TETRAHEDRON, Tet1Co);
215 :
216 : size_t cnt = 0;
217 0 : for(size_t i = 0; i < quadRule.size(); i++, cnt++) {
218 0 : map1.local_to_global(pvPoint[cnt], quadRule.point(i));
219 0 : pvWeight[cnt] = quadRule.weight(i) * map1.sqrt_gram_det(quadRule.point(i));
220 : }
221 :
222 :
223 : DimReferenceMapping<3, 3>& map2 =
224 : ReferenceMappingProvider::get<3,3>(ROID_TETRAHEDRON, Tet2Co);
225 :
226 0 : for(size_t j = 0; j < quadRule.size(); j++, cnt++) {
227 0 : map2.local_to_global(pvPoint[cnt], quadRule.point(j));
228 0 : pvWeight[cnt] = quadRule.weight(j) * map2.sqrt_gram_det(quadRule.point(j));
229 : }
230 0 : m_pvPoint = pvPoint;
231 0 : m_pvWeight = pvWeight;
232 0 : };
233 :
234 0 : GaussQuadraturePyramid::~GaussQuadraturePyramid()
235 : {
236 0 : delete[] m_pvPoint;
237 0 : delete[] m_pvWeight;
238 0 : }
239 :
240 0 : GaussQuadratureOctahedron::GaussQuadratureOctahedron(size_t order)
241 : {
242 0 : GaussQuadratureTetrahedron quadRule = GaussQuadratureTetrahedron(order);
243 :
244 0 : m_order = quadRule.order();
245 0 : m_numPoints = quadRule.size() * 4;
246 0 : position_type* pvPoint = new position_type[m_numPoints];
247 0 : weight_type* pvWeight = new weight_type[m_numPoints];
248 :
249 0 : MathVector<3> Tet1Co[4];
250 : //Tet1Co[0] = MathVector<3>(0,0,0);
251 : //Tet1Co[1] = MathVector<3>(1,1,0);
252 : //Tet1Co[2] = MathVector<3>(0,0,1);
253 : //Tet1Co[3] = MathVector<3>(0,1,0);
254 : Tet1Co[0] = MathVector<3>(0,0,0);
255 : Tet1Co[1] = MathVector<3>(1,0,0);
256 : Tet1Co[2] = MathVector<3>(1,1,0);
257 : Tet1Co[3] = MathVector<3>(0,0,1);
258 :
259 :
260 0 : MathVector<3> Tet2Co[4];
261 : // Tet2Co[0] = MathVector<3>(0,0,0);
262 : // Tet2Co[1] = MathVector<3>(1,0,0);
263 : // Tet2Co[2] = MathVector<3>(1,1,0);
264 : // Tet2Co[3] = MathVector<3>(0,0,1);
265 : Tet2Co[0] = MathVector<3>(0,0,0);
266 : Tet2Co[1] = MathVector<3>(1,1,0);
267 : Tet2Co[2] = MathVector<3>(0,1,0);
268 : Tet2Co[3] = MathVector<3>(0,0,1);
269 :
270 0 : MathVector<3> Tet3Co[4];
271 : // Tet3Co[0] = MathVector<3>(0,0,0);
272 : // Tet3Co[1] = MathVector<3>(1,1,0);
273 : // Tet3Co[2] = MathVector<3>(0,0,-1);
274 : // Tet3Co[3] = MathVector<3>(0,1,0);
275 : Tet3Co[0] = MathVector<3>(0,0,0);
276 : Tet3Co[1] = MathVector<3>(1,1,0);
277 : Tet3Co[2] = MathVector<3>(1,0,0);
278 : Tet3Co[3] = MathVector<3>(0,0,-1);
279 :
280 0 : MathVector<3> Tet4Co[4];
281 : // Tet4Co[0] = MathVector<3>(0,0,0);
282 : // Tet4Co[1] = MathVector<3>(1,0,0);
283 : // Tet4Co[2] = MathVector<3>(1,1,0);
284 : // Tet4Co[3] = MathVector<3>(0,0,-1);
285 : Tet4Co[0] = MathVector<3>(0,0,0);
286 : Tet4Co[1] = MathVector<3>(0,1,0);
287 : Tet4Co[2] = MathVector<3>(1,1,0);
288 : Tet4Co[3] = MathVector<3>(0,0,-1);
289 :
290 :
291 : DimReferenceMapping<3, 3>& map1 =
292 : ReferenceMappingProvider::get<3,3>(ROID_TETRAHEDRON, Tet1Co);
293 :
294 : size_t cnt = 0;
295 0 : for(size_t i = 0; i < quadRule.size(); i++, cnt++) {
296 0 : map1.local_to_global(pvPoint[cnt], quadRule.point(i));
297 0 : pvWeight[cnt] = quadRule.weight(i) * map1.sqrt_gram_det(quadRule.point(i));
298 : }
299 :
300 :
301 : DimReferenceMapping<3, 3>& map2 =
302 : ReferenceMappingProvider::get<3,3>(ROID_TETRAHEDRON, Tet2Co);
303 :
304 0 : for(size_t j = 0; j < quadRule.size(); j++, cnt++) {
305 0 : map2.local_to_global(pvPoint[cnt], quadRule.point(j));
306 0 : pvWeight[cnt] = quadRule.weight(j) * map2.sqrt_gram_det(quadRule.point(j));
307 : }
308 :
309 :
310 : DimReferenceMapping<3, 3>& map3 =
311 : ReferenceMappingProvider::get<3,3>(ROID_TETRAHEDRON, Tet3Co);
312 :
313 0 : for(size_t k = 0; k < quadRule.size(); k++, cnt++) {
314 0 : map3.local_to_global(pvPoint[cnt], quadRule.point(k));
315 0 : pvWeight[cnt] = quadRule.weight(k) * map3.sqrt_gram_det(quadRule.point(k));
316 : }
317 :
318 :
319 : DimReferenceMapping<3, 3>& map4 =
320 : ReferenceMappingProvider::get<3,3>(ROID_TETRAHEDRON, Tet4Co);
321 :
322 0 : for(size_t l = 0; l < quadRule.size(); l++, cnt++) {
323 0 : map4.local_to_global(pvPoint[cnt], quadRule.point(l));
324 0 : pvWeight[cnt] = quadRule.weight(l) * map4.sqrt_gram_det(quadRule.point(l));
325 : }
326 :
327 :
328 0 : m_pvPoint = pvPoint;
329 0 : m_pvWeight = pvWeight;
330 0 : };
331 :
332 0 : GaussQuadratureOctahedron::~GaussQuadratureOctahedron()
333 : {
334 0 : delete[] m_pvPoint;
335 0 : delete[] m_pvWeight;
336 0 : }
337 : }
338 :
|
