1 |
/* $Id$ */ |
2 |
|
3 |
/**************************************************************/ |
4 |
|
5 |
/* assemblage routines: calculate the gradient of nodal data at quadrature points */ |
6 |
|
7 |
/**************************************************************/ |
8 |
|
9 |
/* Copyrights by ACcESS Australia, 2003,2004 */ |
10 |
/* author: gross@access.edu.au */ |
11 |
/* Version: $Id$ */ |
12 |
|
13 |
/**************************************************************/ |
14 |
|
15 |
#include "escript/Data/DataC.h" |
16 |
#include "Common.h" |
17 |
#include "Finley.h" |
18 |
#include "Assemble.h" |
19 |
#include "NodeFile.h" |
20 |
#include "ElementFile.h" |
21 |
#include "Util.h" |
22 |
#ifdef _OPENMP |
23 |
#include <omp.h> |
24 |
#endif |
25 |
/*****************************************************************/ |
26 |
|
27 |
|
28 |
void Finley_Assemble_gradient(Finley_NodeFile* nodes, Finley_ElementFile* elements, |
29 |
escriptDataC* grad_data,escriptDataC* data) { |
30 |
|
31 |
double *local_X=NULL, *local_data=NULL, *dVdv=NULL, *dvdV=NULL, *Vol=NULL, *d_datadv=NULL, *gradS=NULL,*data_array; |
32 |
int numNodes,e,node_offset,*resort_nodes,i,q,type,NS_DOF,NN_DOF,dof_offset; |
33 |
#define NODES 0 |
34 |
#define DOF 1 |
35 |
#define REDUCED_DOF 2 |
36 |
if (nodes==NULL || elements==NULL) return; |
37 |
int NN=elements->ReferenceElement->Type->numNodes; |
38 |
int NS=elements->ReferenceElement->Type->numShapes; |
39 |
int id[NN]; |
40 |
int numDim=nodes->numDim; |
41 |
int data_type=getFunctionSpaceType(data); |
42 |
int numComps=getDataPointSize(data); |
43 |
int numQuad=elements->ReferenceElement->numQuadNodes; |
44 |
for (i=0;i<NN;i++) id[i]=i; |
45 |
|
46 |
/* set some parameter */ |
47 |
|
48 |
if (data_type==FINLEY_NODES) { |
49 |
type=NODES; |
50 |
resort_nodes=id; |
51 |
NN_DOF=elements->ReferenceElement->Type->numNodes; |
52 |
NS_DOF=elements->ReferenceElement->Type->numShapes; |
53 |
gradS=elements->ReferenceElement->dSdv; |
54 |
numNodes=nodes->numNodes; |
55 |
} else if (data_type==FINLEY_DEGREES_OF_FREEDOM) { |
56 |
type=DOF; |
57 |
resort_nodes=id; |
58 |
NN_DOF=elements->ReferenceElement->Type->numNodes; |
59 |
NS_DOF=elements->ReferenceElement->Type->numShapes; |
60 |
gradS=elements->ReferenceElement->dSdv; |
61 |
numNodes=nodes->numDegreesOfFreedom; |
62 |
} else if (data_type==FINLEY_REDUCED_DEGREES_OF_FREEDOM) { |
63 |
type=REDUCED_DOF; |
64 |
resort_nodes=elements->ReferenceElement->Type->linearNodes; |
65 |
NN_DOF=elements->LinearReferenceElement->Type->numNodes; |
66 |
NS_DOF=elements->LinearReferenceElement->Type->numShapes; |
67 |
gradS=elements->LinearReferenceElement->dSdv; |
68 |
numNodes=nodes->reducedNumDegreesOfFreedom; |
69 |
} else { |
70 |
Finley_ErrorCode=TYPE_ERROR; |
71 |
sprintf(Finley_ErrorMsg,"Cannot calculate gradient of data"); |
72 |
} |
73 |
if (getFunctionSpaceType(grad_data)==FINLEY_CONTACT_ELEMENTS_2) { |
74 |
node_offset=NN-NS; |
75 |
dof_offset=NN_DOF-NS_DOF; |
76 |
} else { |
77 |
node_offset=0; |
78 |
dof_offset=0; |
79 |
} |
80 |
|
81 |
/* check the dimensions of interpolated_data and data */ |
82 |
|
83 |
if (numDim!=elements->ReferenceElement->Type->numDim) { |
84 |
Finley_ErrorCode=TYPE_ERROR; |
85 |
sprintf(Finley_ErrorMsg,"Gradient: Spatial and element dimension must match."); |
86 |
} else if (! numSamplesEqual(grad_data,numQuad,elements->numElements)) { |
87 |
Finley_ErrorCode=TYPE_ERROR; |
88 |
sprintf(Finley_ErrorMsg,"illegal number of samples in gradient Data object"); |
89 |
} else if (! numSamplesEqual(data,1,numNodes)) { |
90 |
Finley_ErrorCode=TYPE_ERROR; |
91 |
sprintf(Finley_ErrorMsg,"illegal number of samples of input Data object"); |
92 |
} else if (numDim*numComps!=getDataPointSize(grad_data)) { |
93 |
Finley_ErrorCode=TYPE_ERROR; |
94 |
sprintf(Finley_ErrorMsg,"illegal number of components in gradient data object."); |
95 |
} else if (!isExpanded(grad_data)) { |
96 |
Finley_ErrorCode=TYPE_ERROR; |
97 |
sprintf(Finley_ErrorMsg,"expanded Data object is expected for output data."); |
98 |
} |
99 |
|
100 |
/* now we can start */ |
101 |
|
102 |
if (Finley_ErrorCode==NO_ERROR) { |
103 |
#pragma omp parallel private(local_X,local_data,dvdV,dVdv,Vol,d_datadv) |
104 |
{ |
105 |
local_X=local_data=dVdv=dvdV=Vol=d_datadv=NULL; |
106 |
/* allocation of work arrays */ |
107 |
local_X=THREAD_MEMALLOC(NS*numDim,double); |
108 |
local_data=THREAD_MEMALLOC(NS*numComps,double); |
109 |
dVdv=THREAD_MEMALLOC(numQuad*numDim*numDim,double); |
110 |
dvdV=THREAD_MEMALLOC(numQuad*numDim*numDim,double); |
111 |
Vol=THREAD_MEMALLOC(numQuad,double); |
112 |
d_datadv=THREAD_MEMALLOC(numQuad*numComps*numDim,double); |
113 |
if (!(Finley_checkPtr(local_X) || Finley_checkPtr(dVdv) || Finley_checkPtr(dvdV) || Finley_checkPtr(Vol) || Finley_checkPtr(d_datadv) || Finley_checkPtr(local_data) )) { |
114 |
/* open the element loop */ |
115 |
#pragma omp for private(e,i,q,data_array) schedule(static) |
116 |
for(e=0;e<elements->numElements;e++) { |
117 |
/* gather local coordinates of nodes into local_X: */ |
118 |
Finley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(node_offset,e,NN)]),numDim,nodes->Coordinates,local_X); |
119 |
/* calculate dVdv(i,j,q)=local_X(i,n)*DSDv(n,j,q) */ |
120 |
Finley_Util_SmallMatMult(numDim,numDim*numQuad,dVdv,NS,local_X,elements->ReferenceElement->dSdv); |
121 |
/* dvdV=invert(dVdv) */ |
122 |
Finley_Util_InvertSmallMat(numQuad,numDim,dVdv,dvdV,Vol); |
123 |
/* gather local data into local_data(numComps,NS_DOF): */ |
124 |
switch (type) { |
125 |
case NODES: |
126 |
for (q=0;q<NS_DOF;q++) { |
127 |
i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)]; |
128 |
data_array=getSampleData(data,i); |
129 |
Finley_copyDouble(numComps,data_array,local_data+q*numComps); |
130 |
} |
131 |
break; |
132 |
case DOF: |
133 |
for (q=0;q<NS_DOF;q++) { |
134 |
i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)]; |
135 |
data_array=getSampleData(data,nodes->degreeOfFreedom[i]); |
136 |
Finley_copyDouble(numComps,data_array,local_data+q*numComps); |
137 |
|
138 |
} |
139 |
break; |
140 |
case REDUCED_DOF: |
141 |
for (q=0;q<NS_DOF;q++) { |
142 |
i=elements->Nodes[INDEX2(resort_nodes[dof_offset+q],e,NN)]; |
143 |
data_array=getSampleData(data,nodes->reducedDegreeOfFreedom[i]); |
144 |
Finley_copyDouble(numComps,data_array,local_data+q*numComps); |
145 |
} |
146 |
break; |
147 |
} |
148 |
/* calculate d_datadv(l,i,q)=local_data(l,n)*DSDv(n,i,q) */ |
149 |
Finley_Util_SmallMatMult(numComps,numDim*numQuad,d_datadv,NS_DOF,local_data,gradS); |
150 |
/* calculate grad_data(l,i)=d_datadv(l,i,q)*dvdV(i,k,q) */ |
151 |
Finley_Util_SmallMatSetMult(numQuad,numComps,numDim,getSampleData(grad_data,e),numDim,d_datadv,dvdV); |
152 |
} /* for */ |
153 |
} |
154 |
THREAD_MEMFREE(local_X); |
155 |
THREAD_MEMFREE(dVdv); |
156 |
THREAD_MEMFREE(dvdV); |
157 |
THREAD_MEMFREE(Vol); |
158 |
THREAD_MEMFREE(local_data); |
159 |
THREAD_MEMFREE(d_datadv); |
160 |
} |
161 |
} |
162 |
#undef NODES |
163 |
#undef DOF |
164 |
#undef REDUCED_DOF |
165 |
} |
166 |
/* |
167 |
* $Log$ |
168 |
* Revision 1.4 2004/12/15 07:08:32 jgs |
169 |
* *** empty log message *** |
170 |
* |
171 |
* |
172 |
* |
173 |
*/ |