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__LOCAL_SHAPE_FUNCTION_SET__COMMON__LAGRANGE1D__
34 : #define __H__UG__LIB_DISC__LOCAL_SHAPE_FUNCTION_SET__COMMON__LAGRANGE1D__
35 :
36 : #include "./polynomial1d.h"
37 : #include "common/math/ugmath.h"
38 : #include <vector>
39 :
40 : namespace ug{
41 :
42 : /** Lagrange Polynomial for arbitrary points
43 : *
44 : */
45 : class Lagrange1D
46 : : public Polynomial1D
47 : {
48 : public:
49 : /** constructor for lagrange polynomial i using interpolation points pos
50 : * This constructor creates a lagrange polynomial with interpolation points
51 : * given in pos for the i-th point, i.e. value(pos_i) == 1, value(pos_j) == 0
52 : * for j != i. Therefore, it must hold that 0 <= i < pos.size()
53 : */
54 : Lagrange1D(const size_t i, const std::vector<number>& vPos)
55 : {
56 : // compute coefficients
57 : compute_coeffs(i, vPos);
58 :
59 : // remember positions
60 : m_vPos = vPos;
61 : }
62 :
63 : /// returns the position of the i'th interpolation point
64 : number position(const size_t i) const
65 : {
66 : UG_ASSERT(i < m_vPos.size(), "Invalid index");
67 : return m_vPos[i];
68 : }
69 :
70 : protected:
71 : /// computes the coefficients for passed interpolation points
72 : void compute_coeffs(const size_t i, const std::vector<number>& vPos)
73 : {
74 : // start coefficients
75 : std::vector<number> vStart(1, 1.0);
76 :
77 : // start polynomial
78 : this->set_coefficients(vStart);
79 :
80 : // help polynomial used for each factor
81 : std::vector<number> vFactor(2,1.0);
82 :
83 : // scaling of polynom
84 : number scale = 1.0;
85 :
86 : // fill coefficients
87 : for(size_t j = 0; j < vPos.size(); ++j)
88 : {
89 : if(j == i) continue;
90 :
91 : // set first coefficent to minus position
92 : vFactor[0] = -vPos[j];
93 : // create polynom for (-vPos, 1)
94 : Polynomial1D tmpPol(vFactor);
95 : // multiply
96 : (*this) *= tmpPol;
97 : // multiply scale
98 : scale *= 1./(vPos[i]-vPos[j]);
99 : }
100 :
101 : // multiply by scale
102 : (*this) *= scale;
103 : }
104 :
105 : private:
106 : std::vector<number> m_vPos;
107 : };
108 :
109 : /** EquiDistant Lagrange Function
110 : *
111 : */
112 36 : class EquidistantLagrange1D
113 : : public Polynomial1D
114 : {
115 : public:
116 : /** creates a lagrange polynomial with equidistant interpolation points
117 : * \param[in] i number of interpolation point, where polynom is 1
118 : * \param[in] degree degree of polynom
119 : */
120 36 : EquidistantLagrange1D(const size_t i, const size_t degree)
121 : {
122 : UG_ASSERT(i <= degree, "Only #degree shape functions.");
123 36 : compute_coeffs(i, degree);
124 36 : }
125 :
126 : /// returns the position of the i'th interpolation point
127 : static number position(const size_t i, const size_t degree)
128 : {
129 : UG_ASSERT(i <= degree, "Invalid index");
130 966 : return (number)i/(number)degree;
131 : }
132 :
133 : protected:
134 : /// computes the coefficients for passed interpolation points
135 36 : void compute_coeffs(const int i, const int p)
136 : {
137 : // start coefficients
138 36 : std::vector<number> vStart(1, 1.0);
139 :
140 : // start polynomial
141 36 : this->set_coefficients(vStart);
142 :
143 : // help polynomial used for each factor
144 36 : std::vector<number> vFactor(2, p);
145 :
146 : // scaling of polynom
147 : number scale = 1.0;
148 :
149 : // fill coefficients
150 152 : for(int j = 0; j <= p; ++j)
151 : {
152 116 : if(j == i) continue;
153 :
154 : // set first coefficent to minus position
155 80 : vFactor[0] = -j;
156 : // create polynom for (-j, p)
157 80 : Polynomial1D tmpPol(vFactor);
158 : // multiply
159 80 : (*this) *= tmpPol;
160 : // multiply scale
161 80 : scale *= 1./(i-j);
162 : }
163 :
164 : // multiply by scale
165 : (*this) *= scale;
166 36 : }
167 : };
168 :
169 : /** Truncated EquiDistant Lagrange Function
170 : *
171 : * Creates for given order <tt>p</tt> and iterpolation point <tt>i</tt> the
172 : * polynomial \f[ \prod_{j=0}^{i-1} \frac{x - \frac{j}{p}}{\frac{i}{p} -
173 : * \frac{j}{p}} \f]
174 : */
175 27 : class TruncatedEquidistantLagrange1D
176 : : public Polynomial1D
177 : {
178 : public:
179 : /** creates a lagrange polynomial with equidistant interpolation points
180 : * \param[in] i number of interpolation point, where polynom is 1
181 : * \param[in] degree degree of polynom
182 : */
183 27 : TruncatedEquidistantLagrange1D(const size_t i, const size_t degree)
184 : {
185 : UG_ASSERT(i <= degree, "Only #degree shape functions.");
186 :
187 27 : compute_coeffs(i, degree);
188 27 : }
189 :
190 : /// returns the position of the i'th interpolation point
191 : static number position(const size_t i, const size_t degree)
192 : {
193 : UG_ASSERT(i <= degree, "Invalid index");
194 578 : return (number)i/(number)degree;
195 : }
196 :
197 : protected:
198 : /// computes the coefficients for passed interpolation points
199 27 : void compute_coeffs(const int i, const int p)
200 : {
201 : // start coefficients
202 27 : std::vector<number> vStart(1, 1.0);
203 :
204 : // start polynomial
205 27 : this->set_coefficients(vStart);
206 :
207 : // help polynomial used for each factor
208 27 : std::vector<number> vFactor(2, p);
209 :
210 : // scaling of polynom
211 : number scale = 1.0;
212 :
213 : // fill coefficients
214 57 : for(int j = 0; j < i; ++j)
215 : {
216 : // set first coefficent to minus position
217 30 : vFactor[0] = -j;
218 : // create polynom for (-j, p)
219 30 : Polynomial1D tmpPol(vFactor);
220 : // multiply
221 30 : (*this) *= tmpPol;
222 : // multiply scale
223 30 : scale *= 1./(i-j);
224 : }
225 :
226 : // multiply by scale
227 : (*this) *= scale;
228 27 : }
229 : };
230 :
231 : /** Bounded EquiDistant Lagrange Function
232 : *
233 : * Creates for given order <tt>p</tt>, interpolation point <tt>i</tt> and upper
234 : * bound <tt>0 <= b <= p</tt> the polynomial
235 : * \f[ \prod_{\substack{j=0\\j\neq i}}^{b} \frac{x - \frac{j}{p}}{\frac{i}{p} -
236 : * \frac{j}{p}} \f]
237 : * Thus, it is a polynomial of order b.
238 : */
239 3 : class BoundedEquidistantLagrange1D
240 : : public Polynomial1D
241 : {
242 : public:
243 : /** creates a lagrange polynomial with equidistant interpolation points
244 : * \param[in] i number of interpolation point, where polynom is 1
245 : * \param[in] degree degree of polynom
246 : * \param[in] bound Point until lagrange points are included
247 : */
248 3 : BoundedEquidistantLagrange1D(const size_t i, const size_t degree,
249 : const size_t bound)
250 : {
251 : UG_ASSERT(i <= bound, "Only #bound shape functions.");
252 : UG_ASSERT(bound <= degree, "Only #bound shape functions.");
253 :
254 : // init coefficients
255 3 : compute_coeffs(i, degree, bound);
256 3 : }
257 :
258 : /// returns the position of the i'th interpolation point
259 : static number position(const size_t i, const size_t degree)
260 : {
261 : UG_ASSERT(i <= degree, "Invalid index");
262 : return (number)i/(number)degree;
263 : }
264 :
265 : protected:
266 : /// computes the coefficients for passed interpolation points
267 3 : void compute_coeffs(const int i, const int p, const int b)
268 : {
269 : // start coefficients
270 3 : std::vector<number> vStart(1, 1.0);
271 :
272 : // start polynomial
273 3 : this->set_coefficients(vStart);
274 :
275 : // help polynomial used for each factor
276 3 : std::vector<number> vFactor(2, p);
277 :
278 : // scaling of polynom
279 : number scale = 1.0;
280 :
281 : // fill coefficients
282 8 : for(int j = 0; j <= b; ++j)
283 : {
284 5 : if(j == i) continue;
285 :
286 : // set first coefficent to minus position
287 2 : vFactor[0] = -j;
288 : // create polynom for (-j, p)
289 2 : Polynomial1D tmpPol(vFactor);
290 : // multiply
291 2 : (*this) *= tmpPol;
292 : // multiply scale
293 2 : scale *= 1./(i-j);
294 : }
295 :
296 : // multiply by scale
297 : (*this) *= scale;
298 3 : }
299 : };
300 :
301 :
302 : } // end namespace ug
303 :
304 : #endif /* __H__UG__LIB_DISC__LOCAL_SHAPE_FUNCTION_SET__COMMON__LAGRANGE1D__ */
|
