source: issm/trunk/src/c/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dxt.cpp@ 2549

Last change on this file since 2549 was 2549, checked in by Eric.Larour, 15 years ago

Added multi-threading to InterpFromGridToMeshx.
File size: 5.1 KB
Line 
1/*!\file: InterpFromMeshToMesh2dxt.cpp
2 * \brief "thread" core code for interpolating values from a structured grid.
3 */
4
5#ifdef HAVE_CONFIG_H
6 #include "config.h"
7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
11#include "./InterpFromMeshToMesh2dx.h"
12#include "../shared/shared.h"
13
14#undef __FUNCT__
15#define __FUNCT__ "InterpFromMeshToMesh2dxt"
16void* InterpFromMeshToMesh2dxt(void* vpthread_handle){
17
18 /*gate variables :*/
19 InterpFromMeshToMesh2dxThreadStruct* gate=NULL;
20 pthread_handle* handle=NULL;
21 int my_thread;
22 int num_threads;
23 int i0;
24 int i1;
25 int debug;
26 int nels_data;
27 double* x_data=NULL;
28 double* y_data=NULL;
29 double* index_data=NULL;
30 int nods_prime;
31 double* x_prime=NULL;
32 double* y_prime=NULL;
33 double* data=NULL;
34 int interpolation_type;
35 double default_value;
36 Vec data_prime;
37 double x_prime_min,x_prime_max;
38 double y_prime_min,y_prime_max;
39
40 /*local variables: */
41 int i,j;
42 double x_tria_min,x_tria_max;
43 double y_tria_min,y_tria_max;
44 double area;
45 double area_1,area_2,area_3;
46 double data_value;
47
48 /*recover handle and gate: */
49 handle=(pthread_handle*)vpthread_handle;
50 gate=(InterpFromMeshToMesh2dxThreadStruct*)handle->gate;
51 my_thread=handle->id;
52 num_threads=handle->num;
53
54 /*recover parameters :*/
55 debug= gate->debug;
56 nels_data= gate->nels_data;
57 x_data= gate->x_data;
58 y_data= gate->y_data;
59 index_data= gate->index_data;
60 nods_prime= gate->nods_prime;
61 x_prime= gate->x_prime;
62 y_prime= gate->y_prime;
63 data= gate->data;
64 interpolation_type= gate->interpolation_type;
65 default_value= gate->default_value;
66 data_prime= gate->data_prime;
67 x_prime_min= gate->x_prime_min;
68 x_prime_max= gate->x_prime_max;
69 y_prime_min= gate->y_prime_min;
70 y_prime_max= gate->y_prime_max;
71
72 /*distribute elements across threads :*/
73 PartitionRange(&i0,&i1,nels_data,num_threads,my_thread);
74
75 for (i=i0;i<i1;i++){
76
77 /*display current iteration*/
78 if (debug && (my_thread==0) && fmod((double)(i-i0),(double)100)==0) printf("\b\b\b\b\b\b\b%5.2lf %%",(double)(i-i0)/(i1-i0)*100);
79
80 /*Get extrema coordinates of current elements*/
81 x_tria_min=x_data[(int)index_data[3*i+0]-1]; x_tria_max=x_tria_min;
82 y_tria_min=y_data[(int)index_data[3*i+0]-1]; y_tria_max=y_tria_min;
83
84 for (j=1;j<3;j++){
85 if(x_data[(int)index_data[3*i+j]-1]<x_tria_min) x_tria_min=x_data[(int)index_data[3*i+j]-1];
86 if(x_data[(int)index_data[3*i+j]-1]>x_tria_max) x_tria_max=x_data[(int)index_data[3*i+j]-1];
87 if(y_data[(int)index_data[3*i+j]-1]<y_tria_min) y_tria_min=y_data[(int)index_data[3*i+j]-1];
88 if(y_data[(int)index_data[3*i+j]-1]>y_tria_max) y_tria_max=y_data[(int)index_data[3*i+j]-1];
89 }
90
91 /*if there is no point inside the domain, go to next iteration*/
92 if ( x_prime_max < x_tria_min ) continue;
93 if ( x_prime_min > x_tria_max ) continue;
94 if ( y_prime_max < y_tria_min ) continue;
95 if ( y_prime_min > y_tria_max ) continue;
96
97 /*get area of the current element (Jacobian = 2 * area)*/
98 //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
99 area=x_data[(int)index_data[3*i+1]-1]*y_data[(int)index_data[3*i+2]-1]-y_data[(int)index_data[3*i+1]-1]*x_data[(int)index_data[3*i+2]-1]
100 + x_data[(int)index_data[3*i+0]-1]*y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+0]-1]*x_data[(int)index_data[3*i+1]-1]
101 + x_data[(int)index_data[3*i+2]-1]*y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1]*x_data[(int)index_data[3*i+0]-1];
102
103 /*loop over the prime nodes*/
104 for (j=0;j<nods_prime;j++){
105
106 /*if the current point is not in the triangle, continue*/
107 if ( x_prime[j] < x_tria_min ) continue;
108 if ( x_prime[j] > x_tria_max ) continue;
109 if ( y_prime[j] < y_tria_min ) continue;
110 if ( y_prime[j] > y_tria_max ) continue;
111
112 /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/
113 area_1=((x_prime[j]-x_data[(int)index_data[3*i+2]-1])*(y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+2]-1])
114 - (y_prime[j]-y_data[(int)index_data[3*i+2]-1])*(x_data[(int)index_data[3*i+1]-1]-x_data[(int)index_data[3*i+2]-1]))/area;
115 /*Get second area coordinate =det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/
116 area_2=((x_data[(int)index_data[3*i+0]-1]-x_data[(int)index_data[3*i+2]-1])*(y_prime[j]-y_data[(int)index_data[3*i+2]-1])
117 - (y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1])*(x_prime[j]-x_data[(int)index_data[3*i+2]-1]))/area;
118 /*Get third area coordinate = 1-area1-area2*/
119 area_3=1-area_1-area_2;
120
121 /*is the current point in the current element?*/
122 if (area_1>=0 && area_2>=0 && area_3>=0){
123
124 /*Yes ! compute the value on the point*/
125 if (interpolation_type==1){
126 /*nodal interpolation*/
127 data_value=area_1*data[(int)index_data[3*i+0]-1]+area_2*data[(int)index_data[3*i+1]-1]+area_3*data[(int)index_data[3*i+2]-1];
128 }
129 else{
130 /*element interpolation*/
131 data_value=data[i];
132 }
133 if (isnan(data_value)) data_value=default_value;
134
135 /*insert value and go to the next point*/
136 VecSetValue(data_prime,j,data_value,INSERT_VALUES);
137 }
138 }
139 }
140}
Note: See TracBrowser for help on using the repository browser.