Line data Source code
1 : /*
2 : * Copyright (c) 2010-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__REFERENCE_ELEMENT__REFERENCE_ELEMENT_MAPPING__
34 : #define __H__UG__LIB_DISC__REFERENCE_ELEMENT__REFERENCE_ELEMENT_MAPPING__
35 :
36 : #include <cassert>
37 : #include <iostream>
38 : #include <sstream>
39 : #include "common/common.h"
40 : #include "common/math/ugmath.h"
41 : #include "lib_disc/reference_element/reference_element.h"
42 :
43 : namespace ug{
44 :
45 : extern DebugID DID_REFERENCE_MAPPING;
46 : extern DebugID DID_REFERENCE_MAPPING_GLOB_TO_LOC;
47 :
48 : /**
49 : * This class describes the mapping from a reference element into the real
50 : * (physical) world. The mapping is initialized by the physical positions of
51 : * the vertices of the real world element. The order of those points must be
52 : * given as indicated by the corresponding reference element.
53 : *
54 : * Let \f$R\f$ be the reference element and \f$T\f$ be the element. Then, the
55 : * reference mapping is a mapping:
56 : * \f[
57 : * \phi: R \mapsto T
58 : * \f]
59 : *
60 : * \tparam TRefElem reference element
61 : * \tparam TWorldDim world dimension
62 : */
63 : template <typename TRefElem, int TWorldDim>
64 : class ReferenceMapping
65 : {
66 : public:
67 : /// world dimension (range space dimension)
68 : static const int worldDim = TWorldDim;
69 :
70 : /// reference dimension (domain space dimension)
71 : static const int dim = TRefElem::dim;
72 :
73 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
74 : static const bool isLinear = false;
75 :
76 : public:
77 : /// Default Constructor
78 : ReferenceMapping();
79 :
80 : /// Constructor setting the corners of the element
81 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord);
82 :
83 : /// Constructor setting the corners of the element
84 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord);
85 :
86 : /// refresh mapping for new set of corners
87 : void update(const MathVector<worldDim>* vCornerCoord);
88 :
89 : /// refresh mapping for new set of corners
90 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord);
91 :
92 : /// returns if mapping is affine
93 : bool is_linear() const;
94 :
95 : /// map local coordinate to global coordinate
96 : void local_to_global(MathVector<worldDim>& globPos,
97 : const MathVector<dim>& locPos) const;
98 :
99 : /// map local coordinate to global coordinate for n local positions
100 : void local_to_global(MathVector<worldDim>* vGlobPos,
101 : const MathVector<dim>* vLocPos, size_t n) const;
102 :
103 : /// map local coordinate to global coordinate for a vector of local positions
104 : void local_to_global(std::vector<MathVector<worldDim> >& vGlobPos,
105 : const std::vector<MathVector<dim> >& vLocPos) const;
106 :
107 : /// map global coordinate to local coordinate
108 : void global_to_local(MathVector<dim>& locPos,
109 : const MathVector<worldDim>& globPos,
110 : const size_t maxIter = 1000,
111 : const number tol = 1e-10) const;
112 :
113 : /// map global coordinate to local coordinate for n local positions
114 : void global_to_local(MathVector<dim>* vLocPos,
115 : const MathVector<worldDim>* vGlobPos, size_t n,
116 : const size_t maxIter = 1000,
117 : const number tol = 1e-10) const;
118 :
119 : /// map global coordinate to local coordinate for a vector of local positions
120 : void global_to_local(std::vector<MathVector<dim> >& vLocPos,
121 : const std::vector<MathVector<worldDim> >& vGlobPos,
122 : const size_t maxIter = 1000,
123 : const number tol = 1e-10) const;
124 :
125 : /// returns jacobian
126 : void jacobian(MathMatrix<worldDim, dim>& J,
127 : const MathVector<dim>& locPos) const;
128 :
129 : /// returns jacobian for n local positions
130 : void jacobian(MathMatrix<worldDim, dim>* vJ,
131 : const MathVector<dim>* vLocPos, size_t n) const;
132 :
133 : /// returns jacobian for a vector of local positions
134 : void jacobian(std::vector<MathMatrix<worldDim, dim> >& J,
135 : const std::vector<MathVector<dim> >& vLocPos) const;
136 :
137 : /// returns transposed of jacobian
138 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
139 : const MathVector<dim>& locPos) const;
140 :
141 : /// returns transposed of jacobian for n local positions
142 : void jacobian_transposed(MathMatrix<dim, worldDim>* vJT,
143 : const MathVector<dim>* vLocPos, size_t n) const;
144 :
145 : /// returns transposed of jacobian for a vector of positions
146 : void jacobian_transposed(std::vector<MathMatrix<dim, worldDim> >& vJT,
147 : const std::vector<MathVector<dim> >& vLocPos) const;
148 :
149 : /// returns transposed of the inverse of the jacobian and sqrt of gram determinante
150 : number jacobian_transposed_inverse(MathMatrix<worldDim, dim>& JTInv,
151 : const MathVector<dim>& locPos) const;
152 :
153 : /// returns transposed of the inverse of the jacobian for n local positions
154 : void jacobian_transposed_inverse(MathMatrix<worldDim, dim>* vJTInv,
155 : const MathVector<dim>* vLocPos, size_t n) const;
156 :
157 : /// returns transposed of the inverse of the jacobian for n local positions
158 : void jacobian_transposed_inverse(MathMatrix<worldDim, dim>* vJTInv,
159 : number* vDet,
160 : const MathVector<dim>* vLocPos, size_t n) const;
161 :
162 : /// returns transposed of the inverse of the jacobian for a vector of positions
163 : void jacobian_transposed_inverse(std::vector<MathMatrix<worldDim, dim> >& vJTInv,
164 : const std::vector<MathVector<dim> >& vLocPos) const;
165 :
166 : /// returns transposed of the inverse of the jacobian for a vector of positions
167 : void jacobian_transposed_inverse(std::vector<MathMatrix<worldDim, dim> >& vJTInv,
168 : std::vector<number>& vDet,
169 : const std::vector<MathVector<dim> >& vLocPos) const;
170 :
171 : /// returns the determinate of the jacobian
172 : number sqrt_gram_det(const MathVector<dim>& locPos) const;
173 :
174 : /// returns the determinate of the jacobian for n local positions
175 : void sqrt_gram_det(number* vDet, const MathVector<dim>* vLocPos, size_t n) const;
176 :
177 : /// returns the determinate of the jacobian for a vector of local positions
178 : void sqrt_gram_det(std::vector<number> vDet,
179 : const std::vector<MathVector<dim> >& vLocPos) const;
180 : };
181 :
182 :
183 : ///////////////////////////////////////////////////////////////////////////////
184 : ///////////////////////////////////////////////////////////////////////////////
185 : // Concrete Reference Mappings
186 : ///////////////////////////////////////////////////////////////////////////////
187 : ///////////////////////////////////////////////////////////////////////////////
188 :
189 : /// Base class for Reference mappings helping to implement interface
190 : template <int dim, int worldDim, bool isLinear, typename TImpl>
191 : class BaseReferenceMapping
192 : {
193 : public:
194 : /// returns if mapping is affine
195 : bool is_linear() const {return isLinear;}
196 :
197 : /// map local coordinate to global coordinate for n local positions
198 0 : void local_to_global(MathVector<worldDim>* vGlobPos,
199 : const MathVector<dim>* vLocPos, size_t n) const
200 : {
201 0 : for(size_t ip = 0; ip < n; ++ip)
202 0 : getImpl().local_to_global(vGlobPos[ip], vLocPos[ip]);
203 0 : }
204 :
205 : /// map local coordinate to global coordinate for a vector of local positions
206 0 : void local_to_global(std::vector<MathVector<worldDim> >& vGlobPos,
207 : const std::vector<MathVector<dim> >& vLocPos) const
208 : {
209 : const size_t n = vLocPos.size();
210 0 : vGlobPos.resize(n);
211 0 : local_to_global(&vGlobPos[0], &vLocPos[0], n);
212 0 : }
213 :
214 : /// map global coordinate to local coordinate
215 0 : void global_to_local
216 : (
217 : MathVector<dim>& locPos, ///< (o) for the computed local coordinates (i) specify the initial guess!
218 : const MathVector<worldDim>& globPos, ///< (i) the global coordinates to transform
219 : const size_t maxIter = 1000, ///< (i) maximum number of the Newton iterations
220 : const number tol = 1e-10 ///< (i) tolerance (smalles possible correction in the Newton iterations)
221 : ) const
222 : {
223 : // We apply the Newton's method for the transformation. We assume here that
224 : // the Newton's method converges without the linesearch, and the Jacobian is
225 : // non-singular in the whole iteration process. This is true in particular for
226 : // the bilinear transformation of convex quadrilaterals if the initial guess
227 : // is inside of the element.
228 : MathMatrix<worldDim, dim> J;
229 : MathMatrix<dim, worldDim> JInv;
230 : MathVector<worldDim> dist, compGlobPos, minDist;
231 : MathVector<dim> corr;
232 :
233 : VecScale(minDist, globPos, tol);
234 :
235 0 : for (size_t i = 0; i < maxIter; ++i) {
236 :
237 : // f(x) := \phi(x) - x_{glob}
238 0 : getImpl().local_to_global(compGlobPos, locPos);
239 : VecSubtract(dist, compGlobPos, globPos);
240 :
241 : // Floating-point cancellation protection:
242 0 : if(VecAbsIsLess(dist, minDist))
243 : return;
244 : // REMARK: Note that the tolerance is specified not only to provide the
245 : // sufficient accuracy for the solution but also to protect the iteration
246 : // from the cancellation phenomena in the floating-point arithmetics.
247 : // Computation of the distance involves subtraction which is a reason for
248 : // the loss of precision phenomena. If a small element is located very far
249 : // from the coordinate origin, the accuracy of the local->global transform
250 : // is restricted and this cannot be overcome. After we reach this bound,
251 : // the iterates will oscillate and the defect will not be reduced. We use
252 : // globPos for the reference values.
253 : // Note that this check may not be used as the only stopping criterion:
254 : // globPos may be 0 by specification.
255 :
256 : UG_DLOG(DID_REFERENCE_MAPPING_GLOB_TO_LOC, 2,
257 : "reference_mapping.h: global_to_local() Newton iteration: Iter# "
258 : << i << "; fabs(VecTwoNorm(dist)) = " << fabs(VecTwoNorm(dist)) <<
259 : "; dist = " << dist << "; locPos: " << locPos << std::endl);
260 :
261 : // compute jacobian df/dx = d \phi(x) / dx =: J
262 0 : getImpl().jacobian(J, locPos);
263 :
264 : // solve c -= J^{-1} f
265 0 : LeftInverse(JInv, J);
266 : MatVecMult(corr, JInv, dist);
267 :
268 : // check if tol reached
269 0 : if(VecAbsIsLess(corr, tol))
270 : return;
271 : // REMARK: Note that using the Euclidean or maximum norm directly to dist
272 : // would need tuning of the tolerance for every particular grid: For big
273 : // elements with the diameter of order 1, the tolerance 1e-10 is good accuracy,
274 : // but for a refined grid with the elements of the size of order 1e-5 it
275 : // can be pour. But ||corr|| = ||J^{-1} dist|| is also a norm of dist, and
276 : // it is rescaled according to the element dimensions.
277 : // Besides that, ||corr|| measures whether we can get any further progress
278 : // in the iterations.
279 :
280 : // apply the correction
281 : VecSubtract(locPos, locPos, corr);
282 : }
283 :
284 : // compiler does not know that maxIter will never be 0
285 : // therefore it warns that dist may be uninit'ed;
286 : // as we will throw here anyway, we may as well check that
287 0 : UG_COND_THROW(!maxIter, "Without a single iteration, local-to-global "
288 : "mapping can never converge.");
289 :
290 : UG_DLOG(DID_REFERENCE_MAPPING_GLOB_TO_LOC, 2, "Last JInv:" << std::endl);
291 : for(int i = 0; i < 3; ++i)
292 : {
293 : UG_DLOG(DID_REFERENCE_MAPPING_GLOB_TO_LOC, 2,
294 : JInv(i, 0) << "; " << JInv(i, 1) << "; " << JInv(i, 2) << std::endl);
295 : }
296 :
297 0 : UG_THROW("ReferenceMapping::global_to_local: Newton method did not"
298 : " reach a tolerance "<<tol<<" after "<<maxIter<<" steps. "
299 : "Global Pos: "<<globPos<<", dim: "<<dim<<", worldDim: "<<
300 : worldDim<<", last newton defect: "<<fabs(VecTwoNorm(dist)));
301 : }
302 :
303 : /// map global coordinate to local coordinate for n local positions
304 0 : void global_to_local(MathVector<dim>* vLocPos,
305 : const MathVector<worldDim>* vGlobPos, size_t n,
306 : const size_t maxIter = 1000,
307 : const number tol = 1e-10) const
308 : {
309 : if(isLinear){
310 0 : if(n == 0) return;
311 :
312 : MathMatrix<worldDim, dim> J;
313 : MathMatrix<dim, worldDim> JInv;
314 : MathVector<worldDim> dist, compGlobPos;
315 :
316 : // compute jacobian df/dx = d \phi(x) / dx =: J
317 0 : getImpl().jacobian(J, vLocPos[0]);
318 :
319 : // solve x -= J^{-1} f
320 0 : LeftInverse(JInv, J);
321 :
322 0 : for(size_t ip = 0; ip < n; ++ip)
323 : {
324 : // f(x) := \phi(x) - x_{glob}
325 0 : getImpl().local_to_global(compGlobPos, vLocPos[ip]);
326 0 : VecSubtract(dist, compGlobPos, vGlobPos[ip]);
327 :
328 : MatVecScaleMultAppend(vLocPos[ip], -1.0, JInv, dist);
329 : }
330 : }
331 : else{
332 0 : for(size_t ip = 0; ip < n; ++ip)
333 0 : getImpl().global_to_local(vLocPos[ip], vGlobPos[ip], maxIter, tol);
334 : }
335 : }
336 :
337 : /// map global coordinate to local coordinate for a vector of local positions
338 0 : void global_to_local(std::vector<MathVector<dim> >& vLocPos,
339 : const std::vector<MathVector<worldDim> >& vGlobPos,
340 : const size_t maxIter = 1000,
341 : const number tol = 1e-10) const
342 : {
343 : const size_t n = vGlobPos.size();
344 0 : vLocPos.resize(n);
345 0 : global_to_local(&vLocPos[0], &vGlobPos[0], n, maxIter, tol);
346 0 : }
347 :
348 : /// returns jacobian
349 0 : void jacobian(MathMatrix<worldDim, dim>& J,
350 : const MathVector<dim>& locPos) const
351 : {
352 : MathMatrix<dim, worldDim> JT;
353 0 : getImpl().jacobian_transposed(JT, locPos);
354 : Transpose(J, JT);
355 0 : }
356 :
357 : /// returns jacobian for n local positions
358 0 : void jacobian(MathMatrix<worldDim, dim>* vJ,
359 : const MathVector<dim>* vLocPos, size_t n) const
360 : {
361 : if(isLinear){
362 0 : if(n == 0) return;
363 0 : getImpl().jacobian(vJ[0], vLocPos[0]);
364 0 : for(size_t ip = 1; ip < n; ++ip) vJ[ip] = vJ[0];
365 : }
366 : else {
367 0 : for(size_t ip = 0; ip < n; ++ip)
368 0 : getImpl().jacobian(vJ[ip], vLocPos[ip]);
369 : }
370 : }
371 :
372 : /// returns jacobian for a vector of local positions
373 0 : void jacobian(std::vector<MathMatrix<worldDim, dim> >& vJ,
374 : const std::vector<MathVector<dim> >& vLocPos) const
375 : {
376 : const size_t n = vLocPos.size();
377 0 : vJ.resize(n);
378 0 : jacobian(&vJ[0], &vLocPos[0], n);
379 0 : }
380 :
381 :
382 : /// returns transposed of jacobian for n local positions
383 0 : void jacobian_transposed(MathMatrix<dim, worldDim>* vJT,
384 : const MathVector<dim>* vLocPos, size_t n) const
385 : {
386 : if(isLinear){
387 0 : if(n == 0) return;
388 0 : getImpl().jacobian_transposed(vJT[0], vLocPos[0]);
389 0 : for(size_t ip = 1; ip < n; ++ip) vJT[ip] = vJT[0];
390 : }
391 : else {
392 0 : for(size_t ip = 0; ip < n; ++ip)
393 0 : getImpl().jacobian_transposed(vJT[ip], vLocPos[ip]);
394 : }
395 : }
396 :
397 : /// returns transposed of jacobian for a vector of positions
398 0 : void jacobian_transposed(std::vector<MathMatrix<dim, worldDim> >& vJT,
399 : const std::vector<MathVector<dim> >& vLocPos) const
400 : {
401 : const size_t n = vLocPos.size();
402 0 : vJT.resize(n);
403 0 : jacobian_transposed(&vJT[0], &vLocPos[0], n);
404 0 : }
405 :
406 : /// returns transposed of the inverse of the jacobian
407 0 : number jacobian_transposed_inverse(MathMatrix<worldDim, dim>& JTInv,
408 : const MathVector<dim>& locPos) const
409 : {
410 : MathMatrix<dim, worldDim> JT;
411 0 : getImpl().jacobian_transposed(JT, locPos);
412 : UG_DLOG(DID_REFERENCE_MAPPING, 2, ">>OCT_DISC_DEBUG:: " << "reference_mapping.h: " << "jacobian_transposed_inverse(): JT " << std::endl);
413 : for(int i = 0; i < 3; ++i)
414 : UG_DLOG(DID_REFERENCE_MAPPING, 2, " JT:" << JT(i, 0) << "; " << JT(i, 1) << "; " << JT(i, 2) << std::endl);
415 0 : return RightInverse(JTInv, JT);
416 : }
417 :
418 : /// returns transposed of the inverse of the jacobian for n local positions
419 0 : void jacobian_transposed_inverse(MathMatrix<worldDim, dim>* vJTInv,
420 : number* vDet,
421 : const MathVector<dim>* vLocPos, size_t n) const
422 : {
423 : if(isLinear){
424 0 : if(n == 0) return;
425 0 : vDet[0] = getImpl().jacobian_transposed_inverse(vJTInv[0], vLocPos[0]);
426 0 : for(size_t ip = 1; ip < n; ++ip) vJTInv[ip] = vJTInv[0];
427 0 : for(size_t ip = 1; ip < n; ++ip) vDet[ip] = vDet[0];
428 : }
429 : else {
430 0 : for(size_t ip = 0; ip < n; ++ip)
431 0 : vDet[ip] = getImpl().jacobian_transposed_inverse(vJTInv[ip], vLocPos[ip]);
432 : }
433 : }
434 :
435 : /// returns transposed of the inverse of the jacobian for n local positions
436 0 : void jacobian_transposed_inverse(MathMatrix<worldDim, dim>* vJTInv,
437 : const MathVector<dim>* vLocPos, size_t n) const
438 : {
439 : if(isLinear){
440 0 : if(n == 0) return;
441 0 : getImpl().jacobian_transposed_inverse(vJTInv[0], vLocPos[0]);
442 0 : for(size_t ip = 1; ip < n; ++ip) vJTInv[ip] = vJTInv[0];
443 : }
444 : else {
445 0 : for(size_t ip = 0; ip < n; ++ip)
446 0 : getImpl().jacobian_transposed_inverse(vJTInv[ip], vLocPos[ip]);
447 : }
448 : }
449 :
450 : /// returns transposed of the inverse of the jacobian for a vector of positions
451 0 : void jacobian_transposed_inverse(std::vector<MathMatrix<worldDim, dim> >& vJTInv,
452 : std::vector<number>& vDet,
453 : const std::vector<MathVector<dim> >& vLocPos) const
454 : {
455 : const size_t n = vLocPos.size();
456 0 : vJTInv.resize(n); vDet.resize(n);
457 0 : jacobian_transposed_inverse(&vJTInv[0], &vDet[0], &vLocPos[0], n);
458 0 : }
459 :
460 : /// returns transposed of the inverse of the jacobian for a vector of positions
461 0 : void jacobian_transposed_inverse(std::vector<MathMatrix<worldDim, dim> >& vJTInv,
462 : const std::vector<MathVector<dim> >& vLocPos) const
463 : {
464 : const size_t n = vLocPos.size();
465 0 : vJTInv.resize(n);
466 0 : jacobian_transposed_inverse(&vJTInv[0], &vLocPos[0], n);
467 0 : }
468 :
469 : /// returns the determinate of the jacobian
470 0 : number sqrt_gram_det(const MathVector<dim>& locPos) const
471 : {
472 : MathMatrix<dim, worldDim> JT;
473 0 : getImpl().jacobian_transposed(JT, locPos);
474 0 : return SqrtGramDeterminant(JT);
475 : }
476 :
477 : /// returns the determinate of the jacobian for n local positions
478 0 : void sqrt_gram_det(number* vDet, const MathVector<dim>* vLocPos, size_t n) const
479 : {
480 : if(isLinear){
481 0 : if(n == 0) return;
482 0 : vDet[0] = sqrt_gram_det(vLocPos[0]);
483 0 : for(size_t ip = 1; ip < n; ++ip) vDet[ip] = vDet[0];
484 : }
485 : else {
486 0 : for(size_t ip = 0; ip < n; ++ip)
487 0 : vDet[ip] = sqrt_gram_det(vLocPos[ip]);
488 : }
489 : }
490 :
491 : /// returns the determinate of the jacobian for a vector of local positions
492 0 : void sqrt_gram_det(std::vector<number>& vDet,
493 : const std::vector<MathVector<dim> >& vLocPos) const
494 : {
495 : const size_t n = vLocPos.size();
496 0 : vDet.resize(n);
497 0 : sqrt_gram_det(&vDet[0], &vLocPos[0], n);
498 0 : }
499 :
500 : protected:
501 : /// access to implementation
502 : TImpl& getImpl() {return static_cast<TImpl&>(*this);}
503 :
504 : /// const access to implementation
505 : const TImpl& getImpl() const {return static_cast<const TImpl&>(*this);}
506 : };
507 :
508 : ///////////////////////////////////////////////////////////////////////////////
509 : // Reference Mapping RegularVertex
510 : ///////////////////////////////////////////////////////////////////////////////
511 :
512 : template <int TWorldDim>
513 : class ReferenceMapping<ReferenceVertex, TWorldDim>
514 : : public BaseReferenceMapping<ReferenceVertex::dim, TWorldDim, true,
515 : ReferenceMapping<ReferenceVertex, TWorldDim> >
516 : {
517 : public:
518 : /// world dimension
519 : static const int worldDim = TWorldDim;
520 :
521 : /// reference dimension
522 : static const int dim = ReferenceVertex::dim;
523 :
524 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
525 : static const bool isLinear = true;
526 :
527 : public:
528 : typedef BaseReferenceMapping<ReferenceVertex::dim, TWorldDim, true,
529 : ReferenceMapping<ReferenceVertex, TWorldDim> > base_type;
530 : using base_type::local_to_global;
531 : using base_type::jacobian;
532 : using base_type::jacobian_transposed;
533 : using base_type::jacobian_transposed_inverse;
534 : using base_type::sqrt_gram_det;
535 :
536 : public:
537 : /// Default Constructor
538 : ReferenceMapping() {}
539 :
540 : /// Constructor setting the corners
541 : /// \{
542 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
543 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
544 : /// \}
545 :
546 : /// refresh mapping for new set of corners
547 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
548 : {
549 : UG_ASSERT((int)vCornerCoord.size() >= ReferenceEdge::numCorners,
550 : "ReferenceMapping: to few Corner Coordinates.");
551 : update(&vCornerCoord[0]);
552 : }
553 :
554 : /// update the mapping for a new set of corners
555 : void update(const MathVector<worldDim>* vCornerCoord)
556 : {
557 : co0 = vCornerCoord[0];
558 : }
559 :
560 : /// map local coordinate to global coordinate
561 : void local_to_global(MathVector<worldDim>& globPos,
562 : const MathVector<dim>& locPos) const
563 : {
564 : globPos = co0;
565 : }
566 :
567 : /// returns transposed of jacobian
568 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
569 : const MathVector<dim>& locPos) const
570 : {
571 : //for(int i = 0; i < worldDim; ++i) JT(0,i) = 1;
572 : }
573 :
574 : private:
575 : MathVector<worldDim> co0;
576 : };
577 :
578 : ///////////////////////////////////////////////////////////////////////////////
579 : // Reference Mapping Edge
580 : ///////////////////////////////////////////////////////////////////////////////
581 : template <int TWorldDim>
582 : class ReferenceMapping<ReferenceEdge, TWorldDim>
583 : : public BaseReferenceMapping<ReferenceEdge::dim, TWorldDim, true,
584 : ReferenceMapping<ReferenceEdge, TWorldDim> >
585 : {
586 : public:
587 : /// world dimension
588 : static const int worldDim = TWorldDim;
589 :
590 : /// reference dimension
591 : static const int dim = ReferenceEdge::dim;
592 :
593 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
594 : static const bool isLinear = true;
595 :
596 : public:
597 : typedef BaseReferenceMapping<ReferenceEdge::dim, TWorldDim, true,
598 : ReferenceMapping<ReferenceEdge, TWorldDim> > base_type;
599 : using base_type::local_to_global;
600 : using base_type::jacobian;
601 : using base_type::jacobian_transposed;
602 : using base_type::jacobian_transposed_inverse;
603 : using base_type::sqrt_gram_det;
604 :
605 : public:
606 : /// Default Constructor
607 0 : ReferenceMapping() {}
608 :
609 : /// Constructor setting the corners
610 : /// \{
611 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
612 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
613 : /// \}
614 :
615 : /// refresh mapping for new set of corners
616 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
617 : {
618 : UG_ASSERT((int)vCornerCoord.size() >= ReferenceEdge::numCorners,
619 : "ReferenceMapping: to few Corner Coordinates.");
620 : update(&vCornerCoord[0]);
621 : }
622 :
623 : /// update the mapping for a new set of corners
624 : void update(const MathVector<worldDim>* vCornerCoord)
625 : {
626 : co0 = vCornerCoord[0];
627 : VecSubtract(a10, vCornerCoord[1], vCornerCoord[0]);
628 : }
629 :
630 : /// map local coordinate to global coordinate
631 : void local_to_global(MathVector<worldDim>& globPos,
632 : const MathVector<dim>& locPos) const
633 : {
634 0 : VecScaleAdd(globPos, 1.0, co0, locPos[0], a10);
635 0 : }
636 :
637 : /// returns transposed of jacobian
638 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
639 : const MathVector<dim>& locPos) const
640 : {
641 0 : for(int i = 0; i < worldDim; ++i) JT(0,i) = a10[i];
642 : }
643 :
644 : private:
645 : MathVector<worldDim> co0, a10;
646 : };
647 :
648 : ///////////////////////////////////////////////////////////////////////////////
649 : // Reference Mapping Triangle
650 : ///////////////////////////////////////////////////////////////////////////////
651 : template <int TWorldDim>
652 : class ReferenceMapping<ReferenceTriangle, TWorldDim>
653 : : public BaseReferenceMapping<ReferenceTriangle::dim, TWorldDim, true,
654 : ReferenceMapping<ReferenceTriangle, TWorldDim> >
655 : {
656 : public:
657 : /// world dimension
658 : static const int worldDim = TWorldDim;
659 :
660 : /// reference dimension
661 : static const int dim = ReferenceTriangle::dim;
662 :
663 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
664 : static const bool isLinear = true;
665 :
666 : public:
667 : typedef BaseReferenceMapping<ReferenceTriangle::dim, TWorldDim, true,
668 : ReferenceMapping<ReferenceTriangle, TWorldDim> > base_type;
669 : using base_type::local_to_global;
670 : using base_type::jacobian;
671 : using base_type::jacobian_transposed;
672 : using base_type::jacobian_transposed_inverse;
673 : using base_type::sqrt_gram_det;
674 :
675 : public:
676 : /// Default Constructor
677 0 : ReferenceMapping() {}
678 :
679 : /// Constructor setting the corners
680 : /// \{
681 0 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
682 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
683 : /// \}
684 :
685 : /// refresh mapping for new set of corners
686 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
687 : {
688 : UG_ASSERT((int)vCornerCoord.size() >= ReferenceTriangle::numCorners,
689 : "ReferenceMapping: to few Corner Coordinates.");
690 0 : update(&vCornerCoord[0]);
691 : }
692 :
693 : /// update the mapping for a new set of corners
694 0 : void update(const MathVector<worldDim>* vCornerCoord)
695 : {
696 : co0 = vCornerCoord[0];
697 : VecSubtract(a10, vCornerCoord[1], vCornerCoord[0]);
698 : VecSubtract(a20, vCornerCoord[2], vCornerCoord[0]);
699 0 : }
700 :
701 : /// map local coordinate to global coordinate
702 : void local_to_global(MathVector<worldDim>& globPos,
703 : const MathVector<dim>& locPos) const
704 : {
705 0 : VecScaleAdd(globPos, 1.0, co0, locPos[0], a10, locPos[1], a20);
706 : }
707 :
708 : /// returns transposed of jacobian
709 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
710 : const MathVector<dim>& locPos) const
711 : {
712 0 : for(int i = 0; i < worldDim; ++i) JT(0, i) = a10[i];
713 0 : for(int i = 0; i < worldDim; ++i) JT(1, i) = a20[i];
714 : }
715 :
716 : private:
717 : MathVector<worldDim> co0, a10, a20;
718 :
719 : };
720 : ///////////////////////////////////////////////////////////////////////////////
721 : // Reference Mapping Quadrilateral
722 : ///////////////////////////////////////////////////////////////////////////////
723 :
724 :
725 : template <int TWorldDim>
726 : class ReferenceMapping<ReferenceQuadrilateral, TWorldDim>
727 : : public BaseReferenceMapping<ReferenceQuadrilateral::dim, TWorldDim, false,
728 : ReferenceMapping<ReferenceQuadrilateral, TWorldDim> >
729 : {
730 : public:
731 : /// world dimension
732 : static const int worldDim = TWorldDim;
733 :
734 : /// reference dimension
735 : static const int dim = ReferenceQuadrilateral::dim;
736 :
737 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
738 : static const bool isLinear = false;
739 :
740 : public:
741 : typedef BaseReferenceMapping<ReferenceQuadrilateral::dim, TWorldDim, false,
742 : ReferenceMapping<ReferenceQuadrilateral, TWorldDim> > base_type;
743 : using base_type::local_to_global;
744 : using base_type::jacobian;
745 : using base_type::jacobian_transposed;
746 : using base_type::jacobian_transposed_inverse;
747 : using base_type::sqrt_gram_det;
748 :
749 : public:
750 : /// Default Constructor
751 0 : ReferenceMapping() {}
752 :
753 : /// Constructor setting the corners
754 : /// \{
755 0 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
756 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
757 : /// \}
758 :
759 : /// refresh mapping for new set of corners
760 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
761 : {
762 : UG_ASSERT((int)vCornerCoord.size() >= ReferenceQuadrilateral::numCorners,
763 : "ReferenceMapping: to few Corner Coordinates.");
764 : update(&vCornerCoord[0]);
765 : }
766 :
767 : /// update the mapping for a new set of corners
768 : void update(const MathVector<worldDim>* vCornerCoord)
769 : {
770 0 : for(int co = 0; co < ReferenceQuadrilateral::numCorners; ++co)
771 0 : x[co] = vCornerCoord[co];
772 : }
773 :
774 : /// map local coordinate to global coordinate
775 0 : void local_to_global(MathVector<worldDim>& globPos,
776 : const MathVector<dim>& locPos) const
777 : {
778 0 : const number a = (1.-locPos[0]);
779 0 : const number b = (1.-locPos[1]);
780 :
781 0 : VecScaleAdd(globPos, a*b, x[0],
782 : locPos[0]*b, x[1],
783 : locPos[0]*locPos[1], x[2],
784 : a*locPos[1], x[3]);
785 0 : }
786 :
787 : /// returns transposed of jacobian
788 0 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
789 : const MathVector<dim>& locPos) const
790 : {
791 0 : const number a = 1. - locPos[1];
792 0 : const number b = 1. - locPos[0];
793 :
794 0 : for(int i = 0; i < worldDim; ++i)
795 0 : JT(0, i) = a*(x[1][i] - x[0][i]) + locPos[1]*(x[2][i] - x[3][i]);
796 :
797 0 : for(int i = 0; i < worldDim; ++i)
798 0 : JT(1, i) = b*(x[3][i] - x[0][i]) + locPos[0]*(x[2][i] - x[1][i]);
799 0 : }
800 :
801 : private:
802 : MathVector<worldDim> x[ReferenceQuadrilateral::numCorners];
803 : };
804 :
805 : ///////////////////////////////////////////////////////////////////////////////
806 : // Reference Mapping Tetrahedron
807 : ///////////////////////////////////////////////////////////////////////////////
808 :
809 : template <int TWorldDim>
810 : class ReferenceMapping<ReferenceTetrahedron, TWorldDim>
811 : : public BaseReferenceMapping<ReferenceTetrahedron::dim, TWorldDim, true,
812 : ReferenceMapping<ReferenceTetrahedron, TWorldDim> >
813 : {
814 : public:
815 : /// world dimension
816 : static const int worldDim = TWorldDim;
817 :
818 : /// reference dimension
819 : static const int dim = ReferenceTetrahedron::dim;
820 :
821 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
822 : static const bool isLinear = true;
823 :
824 : public:
825 : typedef BaseReferenceMapping<ReferenceTetrahedron::dim, TWorldDim, true,
826 : ReferenceMapping<ReferenceTetrahedron, TWorldDim> > base_type;
827 : using base_type::local_to_global;
828 : using base_type::jacobian;
829 : using base_type::jacobian_transposed;
830 : using base_type::jacobian_transposed_inverse;
831 : using base_type::sqrt_gram_det;
832 :
833 : public:
834 : /// Default Constructor
835 0 : ReferenceMapping() {}
836 :
837 : /// Constructor setting the corners
838 : /// \{
839 0 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
840 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
841 : /// \}
842 :
843 : /// refresh mapping for new set of corners
844 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
845 : {
846 : UG_ASSERT((int)vCornerCoord.size() >= ReferenceTetrahedron::numCorners,
847 : "ReferenceMapping: to few Corner Coordinates.");
848 0 : update(&vCornerCoord[0]);
849 : }
850 :
851 : /// update the mapping for a new set of corners
852 0 : void update(const MathVector<worldDim>* vCornerCoord)
853 : {
854 : co0 = vCornerCoord[0];
855 : VecSubtract(a10, vCornerCoord[1], vCornerCoord[0]);
856 : VecSubtract(a20, vCornerCoord[2], vCornerCoord[0]);
857 : VecSubtract(a30, vCornerCoord[3], vCornerCoord[0]);
858 0 : }
859 :
860 : /// map local coordinate to global coordinate
861 : void local_to_global(MathVector<worldDim>& globPos,
862 : const MathVector<dim>& locPos) const
863 : {
864 0 : VecScaleAdd(globPos, 1.0, co0,
865 : locPos[0], a10,
866 : locPos[1], a20,
867 : locPos[2], a30);
868 : }
869 :
870 : /// returns transposed of jacobian
871 0 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
872 : const MathVector<dim>& locPos) const
873 : {
874 0 : for(int i = 0; i < worldDim; ++i) JT[0][i] = a10[i];
875 0 : for(int i = 0; i < worldDim; ++i) JT[1][i] = a20[i];
876 0 : for(int i = 0; i < worldDim; ++i) JT[2][i] = a30[i];
877 0 : }
878 :
879 : private:
880 : MathVector<worldDim> co0, a10, a20, a30;
881 : };
882 :
883 : ///////////////////////////////////////////////////////////////////////////////
884 : // Reference Mapping Pyramid
885 : ///////////////////////////////////////////////////////////////////////////////
886 : template <int TWorldDim>
887 : class ReferenceMapping<ReferencePyramid, TWorldDim>
888 : : public BaseReferenceMapping<ReferencePyramid::dim, TWorldDim, false,
889 : ReferenceMapping<ReferencePyramid, TWorldDim> >
890 : {
891 : public:
892 : /// world dimension
893 : static const int worldDim = TWorldDim;
894 :
895 : /// reference dimension
896 : static const int dim = ReferencePyramid::dim;
897 :
898 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
899 : static const bool isLinear = false;
900 :
901 : public:
902 : typedef BaseReferenceMapping<ReferencePyramid::dim, TWorldDim, false,
903 : ReferenceMapping<ReferencePyramid, TWorldDim> > base_type;
904 : using base_type::local_to_global;
905 : using base_type::jacobian;
906 : using base_type::jacobian_transposed;
907 : using base_type::jacobian_transposed_inverse;
908 : using base_type::sqrt_gram_det;
909 :
910 : public:
911 : /// Default Constructor
912 0 : ReferenceMapping() {}
913 :
914 : /// Constructor setting the corners
915 : /// \{
916 0 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
917 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
918 : /// \}
919 :
920 : /// refresh mapping for new set of corners
921 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
922 : {
923 : UG_ASSERT((int)vCornerCoord.size() >= ReferencePyramid::numCorners,
924 : "ReferenceMapping: to few Corner Coordinates.");
925 : update(&vCornerCoord[0]);
926 : }
927 :
928 : /// update the mapping for a new set of corners
929 : void update(const MathVector<worldDim>* vCornerCoord)
930 : {
931 0 : for(int co = 0; co < ReferencePyramid::numCorners; ++co)
932 0 : x[co] = vCornerCoord[co];
933 : }
934 :
935 : /// map local coordinate to global coordinate
936 0 : void local_to_global(MathVector<worldDim>& globPos,
937 : const MathVector<dim>& locPos) const
938 : {
939 0 : const number a = 1.0 - locPos[0];
940 0 : const number b = 1.0 - locPos[1];
941 0 : if (locPos[0] > locPos[1])
942 : {
943 0 : const number a0 = a * b - locPos[2] * b;
944 0 : const number a1 = locPos[0] * b - locPos[2]*locPos[1];
945 0 : const number a2 = locPos[0] * locPos[1] + locPos[2]*locPos[1];
946 0 : const number a3 = a * locPos[1] - locPos[2] * locPos[1];
947 :
948 0 : for(int d = 0; d < worldDim; ++d)
949 0 : globPos[d] = a0*x[0][d]+a1*x[1][d]+a2*x[2][d]
950 0 : +a3*x[3][d]+locPos[2]*x[4][d];
951 : }
952 : else
953 : {
954 0 : const number a0 = a * b - locPos[2] * a;
955 0 : const number a1 = locPos[0] * b - locPos[2]*locPos[0];
956 0 : const number a2 = locPos[0] * locPos[1] + locPos[2]*locPos[0];
957 0 : const number a3 = a * locPos[1] - locPos[2] * locPos[0];
958 0 : for(int d = 0; d < worldDim; ++d)
959 0 : globPos[d] = a0*x[0][d]+a1*x[1][d]+
960 0 : a2*x[2][d]+a3*x[3][d]+locPos[2]*x[4][d];
961 : }
962 0 : }
963 :
964 : /// returns transposed of jacobian
965 0 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
966 : const MathVector<dim>& locPos) const
967 : {
968 : number a[3];
969 0 : for(int d = 0; d < worldDim; ++d)
970 0 : a[d] = x[0][d]-x[1][d]+x[2][d]-x[3][d];
971 :
972 0 : if (locPos[0] > locPos[1])
973 : {
974 0 : for(int d = 0; d < worldDim; ++d)
975 0 : JT(0,d) = x[1][d]-x[0][d]+locPos[1]*a[d];
976 :
977 0 : for(int d = 0; d < worldDim; ++d)
978 0 : JT(1,d) = x[3][d]-x[0][d]+(locPos[0]+locPos[2])*a[d];
979 :
980 0 : for(int d = 0; d < worldDim; ++d)
981 0 : JT(2,d) = x[4][d]-x[0][d]+locPos[1]*a[d];
982 : }
983 : else
984 : {
985 0 : for(int d = 0; d < worldDim; ++d)
986 0 : JT(0,d) = x[1][d]-x[0][d]+(locPos[1]+locPos[2])*a[d];
987 :
988 0 : for(int d = 0; d < worldDim; ++d)
989 0 : JT(1,d) = x[3][d]-x[0][d]+locPos[0]*a[d];
990 :
991 0 : for(int d = 0; d < worldDim; ++d)
992 0 : JT(2,d) = x[4][d]-x[0][d]+locPos[0]*a[d];
993 : }
994 0 : }
995 :
996 : private:
997 : MathVector<worldDim> x[ReferencePyramid::numCorners];
998 : };
999 :
1000 : ///////////////////////////////////////////////////////////////////////////////
1001 : // Reference Mapping Prism
1002 : ///////////////////////////////////////////////////////////////////////////////
1003 :
1004 : template <int TWorldDim>
1005 : class ReferenceMapping<ReferencePrism, TWorldDim>
1006 : : public BaseReferenceMapping<ReferencePrism::dim, TWorldDim, false,
1007 : ReferenceMapping<ReferencePrism, TWorldDim> >
1008 : {
1009 : public:
1010 : /// world dimension
1011 : static const int worldDim = TWorldDim;
1012 :
1013 : /// reference dimension
1014 : static const int dim = ReferencePrism::dim;
1015 :
1016 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
1017 : static const bool isLinear = false;
1018 :
1019 : public:
1020 : typedef BaseReferenceMapping<ReferencePrism::dim, TWorldDim, false,
1021 : ReferenceMapping<ReferencePrism, TWorldDim> > base_type;
1022 : using base_type::local_to_global;
1023 : using base_type::jacobian;
1024 : using base_type::jacobian_transposed;
1025 : using base_type::jacobian_transposed_inverse;
1026 : using base_type::sqrt_gram_det;
1027 :
1028 : public:
1029 : /// Default Constructor
1030 0 : ReferenceMapping() {}
1031 :
1032 : /// Constructor setting the corners
1033 : /// \{
1034 0 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
1035 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
1036 : /// \}
1037 :
1038 : /// refresh mapping for new set of corners
1039 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
1040 : {
1041 : UG_ASSERT((int)vCornerCoord.size() >= ReferencePrism::numCorners,
1042 : "ReferenceMapping: to few Corner Coordinates.");
1043 : update(&vCornerCoord[0]);
1044 : }
1045 :
1046 : /// update the mapping for a new set of corners
1047 : void update(const MathVector<worldDim>* vCornerCoord)
1048 : {
1049 0 : for(int co = 0; co < ReferencePrism::numCorners; ++co)
1050 0 : x[co] = vCornerCoord[co];
1051 : }
1052 :
1053 : /// map local coordinate to global coordinate
1054 0 : void local_to_global(MathVector<worldDim>& globPos,
1055 : const MathVector<dim>& locPos) const
1056 : {
1057 0 : const number a = 1.0 - locPos[0] - locPos[1];
1058 0 : const number b = 1.0 - locPos[2];
1059 0 : const number a0 = a * b;
1060 0 : const number a1 = locPos[0] * b;
1061 0 : const number a2 = locPos[1] * b;
1062 0 : const number a3 = a * locPos[2];
1063 0 : const number a4 = locPos[0] * locPos[2];
1064 0 : const number a5 = locPos[1] * locPos[2];
1065 :
1066 0 : for(int d = 0; d < worldDim; ++d)
1067 0 : globPos[d] = a0*x[0][d]+a1*x[1][d]+a2*x[2][d]+a3*x[3][d]+a4*x[4][d]+a5*x[5][d];
1068 0 : }
1069 :
1070 : /// returns transposed of jacobian
1071 0 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
1072 : const MathVector<dim>& locPos) const
1073 : {
1074 : number a[worldDim], b[worldDim];
1075 :
1076 0 : for(int d = 0; d < worldDim; ++d)
1077 0 : a[d] = x[0][d]-x[1][d]-x[3][d]+x[4][d];
1078 :
1079 0 : for(int d = 0; d < worldDim; ++d)
1080 0 : b[d] = x[0][d]-x[2][d]-x[3][d]+x[5][d];
1081 :
1082 0 : for(int d = 0; d < worldDim; ++d){
1083 0 : JT(0,d) = x[1][d]-x[0][d]+locPos[2]*a[d];
1084 0 : JT(1,d) = x[2][d]-x[0][d]+locPos[2]*b[d];
1085 0 : JT(2,d) = x[3][d]-x[0][d]+locPos[0]*a[d]+locPos[1]*b[d];
1086 : }
1087 0 : }
1088 :
1089 : private:
1090 : MathVector<worldDim> x[ReferencePrism::numCorners];
1091 : };
1092 :
1093 :
1094 : ///////////////////////////////////////////////////////////////////////////////
1095 : // Reference Mapping Hexahedron
1096 : ///////////////////////////////////////////////////////////////////////////////
1097 :
1098 : template <int TWorldDim>
1099 : class ReferenceMapping<ReferenceHexahedron, TWorldDim>
1100 : : public BaseReferenceMapping<ReferenceHexahedron::dim, TWorldDim, false,
1101 : ReferenceMapping<ReferenceHexahedron, TWorldDim> >
1102 : {
1103 : public:
1104 : /// world dimension
1105 : static const int worldDim = TWorldDim;
1106 :
1107 : /// reference dimension
1108 : static const int dim = ReferenceHexahedron::dim;
1109 :
1110 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
1111 : static const bool isLinear = false;
1112 :
1113 : public:
1114 : typedef BaseReferenceMapping<ReferenceHexahedron::dim, TWorldDim, false,
1115 : ReferenceMapping<ReferenceHexahedron, TWorldDim> > base_type;
1116 : using base_type::local_to_global;
1117 : using base_type::jacobian;
1118 : using base_type::jacobian_transposed;
1119 : using base_type::jacobian_transposed_inverse;
1120 : using base_type::sqrt_gram_det;
1121 :
1122 : public:
1123 : /// Default Constructor
1124 0 : ReferenceMapping() {}
1125 :
1126 : /// Constructor setting the corners
1127 : /// \{
1128 0 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
1129 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
1130 : /// \}
1131 :
1132 : /// refresh mapping for new set of corners
1133 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
1134 : {
1135 : UG_ASSERT((int)vCornerCoord.size() >= ReferenceHexahedron::numCorners,
1136 : "ReferenceMapping: to few Corner Coordinates.");
1137 : update(&vCornerCoord[0]);
1138 : }
1139 :
1140 : /// update the mapping for a new set of corners
1141 : void update(const MathVector<worldDim>* vCornerCoord)
1142 : {
1143 0 : for(int co = 0; co < ReferenceHexahedron::numCorners; ++co)
1144 0 : x[co] = vCornerCoord[co];
1145 : }
1146 :
1147 : /// map local coordinate to global coordinate
1148 0 : void local_to_global(MathVector<worldDim>& globPos,
1149 : const MathVector<dim>& locPos) const
1150 : {
1151 0 : const number a = 1.0 - locPos[0];
1152 0 : const number b = 1.0 - locPos[1];
1153 0 : const number c = 1.0 - locPos[2];
1154 0 : const number a0 = a * b * c;
1155 0 : const number a1 = locPos[0] * b * c;
1156 0 : const number a2 = locPos[0] * locPos[1] * c;
1157 0 : const number a3 = a * locPos[1] * c;
1158 0 : const number a4 = a * b * locPos[2];
1159 0 : const number a5 = locPos[0] * b * locPos[2];
1160 0 : const number a6 = locPos[0] * locPos[1] * locPos[2];
1161 0 : const number a7 = a * locPos[1] * locPos[2];
1162 :
1163 0 : for(int d = 0; d < worldDim; ++d){
1164 0 : globPos[d] = a0*x[0][d]+a1*x[1][d]+a2*x[2][d]+a3*x[3][d]+
1165 0 : a4*x[4][d]+a5*x[5][d]+a6*x[6][d]+a7*x[7][d];
1166 : }
1167 0 : }
1168 :
1169 : /// returns transposed of jacobian
1170 0 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
1171 : const MathVector<dim>& locPos) const
1172 : {
1173 : number a0,a1,a2,a3;
1174 0 : const number a = 1.0 - locPos[0];
1175 0 : const number b = 1.0 - locPos[1];
1176 0 : const number c = 1.0 - locPos[2];
1177 0 : a0 = b * c;
1178 0 : a1 = locPos[1] * c;
1179 0 : a2 = locPos[1] * locPos[2];
1180 0 : a3 = b * locPos[2];
1181 0 : for(int d = 0; d < worldDim; ++d)
1182 0 : JT(0,d) = a0*(x[1][d]-x[0][d])+a1*(x[2][d]-x[3][d])
1183 0 : + a2*(x[6][d]-x[7][d])+a3*(x[5][d]-x[4][d]);
1184 :
1185 0 : a0 = a * c;
1186 0 : a1 = locPos[0] * c;
1187 0 : a2 = locPos[0] * locPos[2];
1188 0 : a3 = a * locPos[2];
1189 0 : for(int d = 0; d < worldDim; ++d)
1190 0 : JT(1,d) = a0*(x[3][d]-x[0][d])+a1*(x[2][d]-x[1][d])
1191 0 : + a2*(x[6][d]-x[5][d])+a3*(x[7][d]-x[4][d]);
1192 :
1193 0 : a0 = a * b;
1194 0 : a1 = locPos[0] * b;
1195 0 : a2 = locPos[0] * locPos[1];
1196 0 : a3 = a * locPos[1];
1197 0 : for(int d = 0; d < worldDim; ++d)
1198 0 : JT(2,d) = a0*(x[4][d]-x[0][d])+a1*(x[5][d]-x[1][d])
1199 0 : + a2*(x[6][d]-x[2][d])+a3*(x[7][d]-x[3][d]);
1200 0 : }
1201 :
1202 : private:
1203 : MathVector<worldDim> x[ReferenceHexahedron::numCorners];
1204 : };
1205 :
1206 : ///////////////////////////////////////////////////////////////////////////////
1207 : // Reference Mapping Octahedron
1208 : ///////////////////////////////////////////////////////////////////////////////
1209 : template <int TWorldDim>
1210 : class ReferenceMapping<ReferenceOctahedron, TWorldDim>
1211 : : public BaseReferenceMapping<ReferenceOctahedron::dim, TWorldDim, false,
1212 : ReferenceMapping<ReferenceOctahedron, TWorldDim> >
1213 : {
1214 : public:
1215 : /// world dimension
1216 : static const int worldDim = TWorldDim;
1217 :
1218 : /// reference dimension
1219 : static const int dim = ReferenceOctahedron::dim;
1220 :
1221 : /// flag if mapping is linear (i.e. Jacobian does not depend on x)
1222 : static const bool isLinear = false;
1223 :
1224 : public:
1225 : typedef BaseReferenceMapping<ReferenceOctahedron::dim, TWorldDim, false,
1226 : ReferenceMapping<ReferenceOctahedron, TWorldDim> > base_type;
1227 : using base_type::local_to_global;
1228 : using base_type::jacobian;
1229 : using base_type::jacobian_transposed;
1230 : using base_type::jacobian_transposed_inverse;
1231 : using base_type::sqrt_gram_det;
1232 :
1233 : public:
1234 : /// Default Constructor
1235 0 : ReferenceMapping() {}
1236 :
1237 : /// Constructor setting the corners
1238 : /// \{
1239 0 : ReferenceMapping(const MathVector<worldDim>* vCornerCoord) {update(vCornerCoord);}
1240 : ReferenceMapping(const std::vector<MathVector<worldDim> >& vCornerCoord) {update(vCornerCoord);}
1241 : /// \}
1242 :
1243 : /// refresh mapping for new set of corners
1244 : void update(const std::vector<MathVector<worldDim> >& vCornerCoord)
1245 : {
1246 : UG_ASSERT((int)vCornerCoord.size() >= ReferenceOctahedron::numCorners,
1247 : "ReferenceMapping: to few Corner Coordinates.");
1248 : update(&vCornerCoord[0]);
1249 : }
1250 :
1251 : /// update the mapping for a new set of corners
1252 : void update(const MathVector<worldDim>* vCornerCoord)
1253 : {
1254 0 : for(int co = 0; co < ReferenceOctahedron::numCorners; ++co)
1255 0 : x[co] = vCornerCoord[co];
1256 : }
1257 :
1258 : /// map local coordinate to global coordinate
1259 0 : void local_to_global(MathVector<worldDim>& globPos,
1260 : const MathVector<dim>& locPos) const
1261 : {
1262 : // the octahedral shape functions correspond to the tetrahedral ones,
1263 : // but locally piecewise linear on a subdivision of the octahedron
1264 : // into 4 sub-tetrahedrons.
1265 0 : const number z_sgn = (locPos[2] < 0) ? -1.0 : 1.0;
1266 0 : const number a0 = (locPos[2] < 0) ? -locPos[2] : 0.0;
1267 0 : const number a5 = (locPos[2] < 0) ? 0.0 : locPos[2];
1268 :
1269 0 : if (locPos[0] > locPos[1])
1270 : {
1271 0 : const number a1 = 1.0 - locPos[0] - z_sgn * locPos[2];
1272 0 : const number a2 = locPos[0] - locPos[1];
1273 : const number a3 = locPos[1];
1274 : const number a4 = 0.0;
1275 :
1276 0 : for(int d = 0; d < worldDim; ++d)
1277 0 : globPos[d] = a0*x[0][d]+a1*x[1][d]+a2*x[2][d]
1278 0 : +a3*x[3][d]+a4*x[4][d]+a5*x[5][d];
1279 : }
1280 : else
1281 : {
1282 0 : const number a1 = 1.0 - locPos[1] - z_sgn * locPos[2];
1283 : const number a2 = 0.0;
1284 : const number a3 = locPos[0];
1285 0 : const number a4 = locPos[1] - locPos[0];
1286 :
1287 0 : for(int d = 0; d < worldDim; ++d)
1288 0 : globPos[d] = a0*x[0][d]+a1*x[1][d]+a2*x[2][d]
1289 0 : +a3*x[3][d]+a4*x[4][d]+a5*x[5][d];
1290 : }
1291 0 : }
1292 :
1293 : /// returns transposed of jacobian
1294 0 : void jacobian_transposed(MathMatrix<dim, worldDim>& JT,
1295 : const MathVector<dim>& locPos) const
1296 : {
1297 : // the octahedral shape functions correspond to the tetrahedral ones,
1298 : // but locally piecewise linear on a subdivision of the octahedron
1299 : // into 4 sub-tetrahedrons.
1300 0 : const number z_sgn = (locPos[2] < 0) ? -1.0 : 1.0;
1301 0 : const number Da0 = (locPos[2] < 0) ? -1.0 : 0.0;
1302 0 : const number Da5 = (locPos[2] < 0) ? 0.0 : 1.0;
1303 :
1304 0 : if (locPos[0] > locPos[1])
1305 : {
1306 0 : for(int d = 0; d < worldDim; ++d)
1307 0 : JT(0,d) = -x[1][d]+x[2][d];
1308 :
1309 0 : for(int d = 0; d < worldDim; ++d)
1310 0 : JT(1,d) = -x[2][d]+x[3][d];
1311 :
1312 0 : for(int d = 0; d < worldDim; ++d)
1313 0 : JT(2,d) = Da0*x[0][d] - z_sgn*x[1][d] + Da5*x[5][d];
1314 : }
1315 : else
1316 : {
1317 0 : for(int d = 0; d < worldDim; ++d)
1318 0 : JT(0,d) = x[3][d]-x[4][d];
1319 :
1320 0 : for(int d = 0; d < worldDim; ++d)
1321 0 : JT(1,d) = -x[1][d]+x[4][d];
1322 :
1323 0 : for(int d = 0; d < worldDim; ++d)
1324 0 : JT(2,d) = Da0*x[0][d] - z_sgn*x[1][d] + Da5*x[5][d];
1325 : }
1326 0 : }
1327 :
1328 : private:
1329 : MathVector<worldDim> x[ReferenceOctahedron::numCorners];
1330 : };
1331 :
1332 :
1333 : } // end namespace ug
1334 :
1335 : #endif /* __H__UG__LIB_DISC__REFERENCE_ELEMENT__REFERENCE_ELEMENT_MAPPING__ */
|
