【讨论】域内点插值

hhf_cas

我已经有一系列平面单元（1255个节点，非规则，含三节点三角单元、四节点四边形单元）及节点的计算结果，需将这些节点计算结果插值到一系列规则正方形单元的的节点上，请问如何实现？力求计算量较小。
补充：我这后续软件只能识别方格子上的信息，边界可以忽略

zccbest

还不如重新分网,重新计算

hhf_cas

2# zccbest 对于一个不是方正的结构（比如说圆），根本难以画出规则格子。我这后续软件只能识别方格子上的信息，边界可以忽略

zccbest

这个比较复杂吧,先要确定后来的节点坐落在原先哪个单元内,根据原来单元的插值函数求解单元内任意一点的函数值

hhf_cas

3# hhf_cas 确实是需要判断，所以我想求教达人，有无比较快速的算法？

hhf_cas

4# zccbest 确实是需要判断，所以我想求教达人，有无比较快速的算法

chentao807

方形的格子你怎么画？以我的理解，你的几何很难划分成结构单元。

hhf_cas

7# chentao807 自己用程序生成格子，落在域内的点插值，域外的赋0

fanxiaochen

这个要用到搜索算法，计算量应该少不了，同时还要保证精度，建议参考无网格用的插值方法

hhf_cas

9# fanxiaochen 我又这个想法，只是计算速度有待实践。估计边界上的精度不会太好

pasuka

我已经有一系列平面单元（1255个节点，非规则，含三节点三角单元、四节点四边形单元）及节点的计算结果，需将这些节点计算结果插值到一系列规则正方形单元的的节点上，请问如何实现？力求计算量较小。
补充：我这后 ...
hhf_cas 发表于 2011-7-12 09:38

既然有这个要求，为何不用边界元或者无网格方法呢？
甚至有限差分法都可以，如果边界不复杂

caoer

your idea is doable.

here is my interpolation class written in C++.
No public variable exists, so the header file is ignored.

1. #include "polynomial_interpolants.hpp"
   
2. 
   
3. #include "../general/spectral_common.hpp"
   
4. 
   
5. 
   
6. 
   
7. //algorithem 30
   
8. 
   
9. void PolynomialInterpolants::BraycentricWeigths(const int N,
   
10. 
    
11.                                                 const double *x_vec,
    
12. 
    
13.                                                 double  *w_vec)
    
14. 
    
15. {
    
16. 
    
17.     int j,k;
    
18. 
    
19. 
    
20. 
    
21.     for (j = 0; j<=N; j++)
    
22. 
    
23.         *(w_vec+j) = 1.0;
    
24. 
    
25. 
    
26. 
    
27.     for (j=1; j<=N; j++)
    
28. 
    
29.     {
    
30. 
    
31.         for (k = 0; k<j; k++)
    
32. 
    
33.         {
    
34. 
    
35.             *(w_vec+k) = *(w_vec+k) * (*(x_vec+k) - *(x_vec+j));
    
36. 
    
37.             *(w_vec+j) = *(w_vec+j) * (*(x_vec+j) - *(x_vec+k));
    
38. 
    
39.         }
    
40. 
    
41.     }
    
42. 
    
43.     for (j=0; j<=N; j++)
    
44. 
    
45.         *(w_vec+j) = 1/ *(w_vec+j);
    
46. 
    
47. };
    
48. 
    
49. 
    
50. 
    
51. 
    
52. 
    
53. 
    
54. 
    
55. // Algorithm 31
    
56. 
    
57. double PolynomialInterpolants::LagrangeInterpolation(const unsigned int N,
    
58. 
    
59.                                                      const double x,
    
60. 
    
61.                                                      const double *x_vec,
    
62. 
    
63.                                                      const double *f_vec,
    
64. 
    
65.                                                      const double *w_vec)
    
66. 
    
67. {
    
68. 
    
69.     double t;
    
70. 
    
71.     double numerator = 0.0, denominator = 0.0;
    
72. 
    
73.     SpectralCommon *pSC = new SpectralCommon();
    
74. 
    
75. 
    
76. 
    
77. //    std::cout<<"N = "<< N << std::endl;
    
78. 
    
79. //    for (unsigned int i=0; i<= N; i++)
    
80. 
    
81. //        std::cout<<x_vec[i] <<" , "<< f_vec[i]<<" , "<< w_vec[i]<<std::endl;
    
82. 
    
83. 
    
84. 
    
85.     for (unsigned int j = 0; j<=N; j++)
    
86. 
    
87.     {
    
88. 
    
89.         if (pSC->AlmostEqual(x, x_vec[j]))
    
90. 
    
91.         {
    
92. 
    
93.             return f_vec[j];
    
94. 
    
95.         }
    
96. 
    
97.         else
    
98. 
    
99.         {
    
100. 
     
101.             t = w_vec[j]/(x - x_vec[j]);
     
102. 
     
103.             numerator   = numerator + t* f_vec[j];
     
104. 
     
105.             denominator = denominator + t;
     
106. 
     
107.         }
     
108. 
     
109.     }
     
110. 
     
111. 
     
112. 
     
113.     delete pSC;
     
114. 
     
115.     return numerator/denominator;
     
116. 
     
117. }
     
118. 
     
119. 
     
120. 
     
121. 
     
122. 
     
123. 
     
124. 
     
125. // Algorithm 32
     
126. 
     
127. void PolynomialInterpolants::PolynomialInterpolationMatrix(const int N,
     
128. 
     
129.                                                            const int M,
     
130. 
     
131.                                                            const double *x_vec,
     
132. 
     
133.                                                            const double *w_vec,
     
134. 
     
135.                                                            const double *xi_vec,
     
136. 
     
137.                                                            double *T_mtx)
     
138. 
     
139. {
     
140. 
     
141.     bool rowHasMatch;
     
142. 
     
143.     int k, j;
     
144. 
     
145.     double s, t;
     
146. 
     
147.     SpectralCommon *pSC = new SpectralCommon();
     
148. 
     
149. 
     
150. 
     
151.     for (k=0; k<=M; k++)
     
152. 
     
153.     {
     
154. 
     
155.         rowHasMatch = false;
     
156. 
     
157.         for (j=0; j<=N; j++)
     
158. 
     
159.         {
     
160. 
     
161.             *(T_mtx+k*(N+1) + j)=0;
     
162. 
     
163. 
     
164. 
     
165.             if (pSC->AlmostEqual(*(xi_vec+k), *(x_vec+j)))
     
166. 
     
167.             {
     
168. 
     
169.                 rowHasMatch = true;
     
170. 
     
171.                 *(T_mtx+k*(N+1)+j) = 1;
     
172. 
     
173.             } //if
     
174. 
     
175.         } // for j
     
176. 
     
177.         if (!rowHasMatch)
     
178. 
     
179.         {
     
180. 
     
181.             s = 0;
     
182. 
     
183.             for(j=0; j<=N; j++)
     
184. 
     
185.             {
     
186. 
     
187.                 t = *(w_vec+j)/(*(xi_vec+k)- *(x_vec+j));
     
188. 
     
189.                 *(T_mtx+k*(N+1)+j) = t;
     
190. 
     
191.                 s = s + t;
     
192. 
     
193.             }
     
194. 
     
195.             for (j = 0; j<=N; j++)
     
196. 
     
197.                 *(T_mtx +k*(N+1)+j) = *(T_mtx +k*(N+1)+j) / s;
     
198. 
     
199. 
     
200. 
     
201.         } // if
     
202. 
     
203.     } //for k
     
204. 
     
205.     delete pSC;
     
206. 
     
207. }
     
208. 
     
209. 
     
210. 
     
211. 
     
212. 
     
213. 
     
214. 
     
215. // Algorithm 33
     
216. 
     
217. void PolynomialInterpolants::InterpolateToNewpoints(const int M,
     
218. 
     
219.                                                     const int N,
     
220. 
     
221.                                                     const double *T_mtx,
     
222. 
     
223.                                                     const double *f_vec,
     
224. 
     
225.                                                     double *fInterp_vec)
     
226. 
     
227. {
     
228. 
     
229.     int i,j;
     
230. 
     
231.     double t;
     
232. 
     
233.     for(i =0; i<=M; i++)
     
234. 
     
235.     {
     
236. 
     
237.         t = 0;
     
238. 
     
239.         for(j=0; j<=N; j++)
     
240. 
     
241.             t = t+ *(T_mtx +i*(N+1)+j) * *(f_vec +j);
     
242. 
     
243.         *(fInterp_vec+i) = t;
     
244. 
     
245.     }
     
246. 
     
247. }
     
248. 
     
249. 
     
250. 
     
251. 
     
252. 
     
253. // Algorithm 34
     
254. 
     
255. void PolynomialInterpolants::LagrangeInterpolatingPolynomials(const int N,
     
256. 
     
257.                                                               const double x,
     
258. 
     
259.                                                               const double *x_vec,
     
260. 
     
261.                                                               const double *w_vec,
     
262. 
     
263.                                                               double *l_vec)
     
264. 
     
265. {
     
266. 
     
267.     bool xMatchesNode= false;
     
268. 
     
269.     double s, t;
     
270. 
     
271.     int j;
     
272. 
     
273.     SpectralCommon *pSC = new SpectralCommon();
     
274. 
     
275. 
     
276. 
     
277.     for (j =0; j<=N; j++)
     
278. 
     
279.     {
     
280. 
     
281.         *(l_vec+j) = 0.0;
     
282. 
     
283.         if (pSC->AlmostEqual(x, *(x_vec+j)))
     
284. 
     
285.         {
     
286. 
     
287.             *(l_vec+j) = 1.0;
     
288. 
     
289.             xMatchesNode = true;
     
290. 
     
291.         }
     
292. 
     
293.     } // for j
     
294. 
     
295. 
     
296. 
     
297.     if (!xMatchesNode)
     
298. 
     
299.     {
     
300. 
     
301.         s = 0.0;
     
302. 
     
303.         for (j = 0; j <= N; j++)
     
304. 
     
305.         {
     
306. 
     
307.             t = *(w_vec+j)/(x - *(x_vec+j));
     
308. 
     
309.             *(l_vec+j) = t;
     
310. 
     
311.             s = s + t;
     
312. 
     
313.         }
     
314. 
     
315.         for(j = 0; j<= N; j++)
     
316. 
     
317.             *(l_vec+j) = *(l_vec+j) / s;
     
318. 
     
319.     } //if
     
320. 
     
321. 
     
322. 
     
323.     delete pSC;
     
324. 
     
325. }

复制代码

hhf_cas

12# caoer Thank you so much. I'm deeply moved by your enthusiasm.