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

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

New multi-threading capability enabled for InterpFromMeshToMesh2d module.
File size: 5.3 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 int step;
48
49 /*recover handle and gate: */
50 handle=(pthread_handle*)vpthread_handle;
51 gate=(InterpFromMeshToMesh2dxThreadStruct*)handle->gate;
52 my_thread=handle->id;
53 num_threads=handle->num;
54
55 /*recover parameters :*/
56 debug= gate->debug;
57 nels_data= gate->nels_data;
58 x_data= gate->x_data;
59 y_data= gate->y_data;
60 index_data= gate->index_data;
61 nods_prime= gate->nods_prime;
62 x_prime= gate->x_prime;
63 y_prime= gate->y_prime;
64 data= gate->data;
65 interpolation_type= gate->interpolation_type;
66 default_value= gate->default_value;
67 data_prime= gate->data_prime;
68 x_prime_min= gate->x_prime_min;
69 x_prime_max= gate->x_prime_max;
70 y_prime_min= gate->y_prime_min;
71 y_prime_max= gate->y_prime_max;
72
73 /*distribute elements across threads :*/
74 step=(int)floor((double)nels_data/(double)num_threads);
75 for(i=0;i<(my_thread+1);i++){
76 i0=i*step;
77 if(i==(num_threads-1))i1=nels_data;
78 else i1=i0+step;
79 }
80
81 for (i=i0;i<i1;i++){
82
83 /*display current iteration*/
84 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);
85
86 /*Get extrema coordinates of current elements*/
87 x_tria_min=x_data[(int)index_data[3*i+0]-1]; x_tria_max=x_tria_min;
88 y_tria_min=y_data[(int)index_data[3*i+0]-1]; y_tria_max=y_tria_min;
89
90 for (j=1;j<3;j++){
91 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];
92 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];
93 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];
94 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];
95 }
96
97 /*if there is no point inside the domain, go to next iteration*/
98 if ( x_prime_max < x_tria_min ) continue;
99 if ( x_prime_min > x_tria_max ) continue;
100 if ( y_prime_max < y_tria_min ) continue;
101 if ( y_prime_min > y_tria_max ) continue;
102
103 /*get area of the current element (Jacobian = 2 * area)*/
104 //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
105 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]
106 + 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]
107 + 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];
108
109 /*loop over the prime nodes*/
110 for (j=0;j<nods_prime;j++){
111
112 /*if the current point is not in the triangle, continue*/
113 if ( x_prime[j] < x_tria_min ) continue;
114 if ( x_prime[j] > x_tria_max ) continue;
115 if ( y_prime[j] < y_tria_min ) continue;
116 if ( y_prime[j] > y_tria_max ) continue;
117
118 /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/
119 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])
120 - (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;
121 /*Get second area coordinate =det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/
122 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])
123 - (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;
124 /*Get third area coordinate = 1-area1-area2*/
125 area_3=1-area_1-area_2;
126
127 /*is the current point in the current element?*/
128 if (area_1>=0 && area_2>=0 && area_3>=0){
129
130 /*Yes ! compute the value on the point*/
131 if (interpolation_type==1){
132 /*nodal interpolation*/
133 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];
134 }
135 else{
136 /*element interpolation*/
137 data_value=data[i];
138 }
139 if (isnan(data_value)) data_value=default_value;
140
141 /*insert value and go to the next point*/
142 VecSetValue(data_prime,j,data_value,INSERT_VALUES);
143 }
144 }
145 }
146}
Note: See TracBrowser for help on using the repository browser.