/[escript]/trunk/esys2/finley/src/finleyC/Assemble_RobinCondition.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/Assemble_RobinCondition.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 11 months ago) by jgs
File MIME type: text/plain
File size: 10327 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* assembles a system of numEq natural boundary condition into the stiffness matrix S and right hand side F: */
6
7 /* d*u = y */
8
9 /* u has numComp components. */
10
11 /* Shape of the coefficients: */
12
13 /* d = numEqu x numComp */
14 /* y = numEqu */
15
16 /* The coefficients d and y have to be defined on the integartion points on face elements or not present (=NULL). */
17
18 /* S and F have to be initialized before the routine is called. S or F can be NULL. In this case the left or */
19 /* the right hand side of the natural boundary conditions is not processed. */
20
21
22 /**************************************************************/
23
24 /* Copyrights by ACcESS Australia, 2003,2004 */
25 /* author: gross@access.edu.au */
26 /* Version: $Id$ */
27
28 /**************************************************************/
29
30 #include "escript/Data/DataC.h"
31 #include "Assemble.h"
32 #include "NodeFile.h"
33 #include "ElementFile.h"
34 #include "Finley.h"
35 #include "Util.h"
36 #ifdef _OPENMP
37 #include <omp.h>
38 #endif
39
40 /**************************************************************/
41
42 void Finley_Assemble_RobinCondition(Finley_NodeFile* nodes,Finley_ElementFile* elements,Finley_SystemMatrix* S, escriptDataC* F, escriptDataC* d, escriptDataC* y, Finley_Assemble_handelShapeMissMatch handelShapeMissMatchForEM) {
43 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Area=NULL;
44 double time0;
45 Assemble_Parameters p;
46 maybelong *index_row=NULL,*index_col=NULL;
47 int dimensions[ESCRIPT_MAX_DATA_RANK],e,q,color;
48
49 if (nodes==NULL || elements==NULL) return;
50 if (S==NULL && isEmpty(F)) return;
51
52 /* set all parameters in p*/
53 Assemble_getAssembleParameters(nodes,elements,S,F,&p);
54 if (Finley_ErrorCode!=NO_ERROR) return;
55
56 /* get a functionspace */
57 int funcspace=UNKNOWN;
58 updateFunctionSpaceType(funcspace,d);
59 updateFunctionSpaceType(funcspace,y);
60 if (funcspace==UNKNOWN) return; /* all data are empty */
61
62 /* check if all function spaces are the same */
63
64 if (! functionSpaceTypeEqual(funcspace,d) ) {
65 Finley_ErrorCode=TYPE_ERROR;
66 sprintf(Finley_ErrorMsg,"unexpected function space typ for coefficient d");
67 }
68 if (! functionSpaceTypeEqual(funcspace,y) ) {
69 Finley_ErrorCode=TYPE_ERROR;
70 sprintf(Finley_ErrorMsg,"unexpected function space typ for coefficient y");
71 }
72
73 /* check if all function spaces are the same */
74
75 if (! numSamplesEqual(d,p.numQuad,elements->numElements) ) {
76 Finley_ErrorCode=TYPE_ERROR;
77 sprintf(Finley_ErrorMsg,"sample points of coefficient d don't match (%d,%d)",p.numQuad,elements->numElements);
78 }
79
80 if (! numSamplesEqual(y,p.numQuad,elements->numElements) ) {
81 Finley_ErrorCode=TYPE_ERROR;
82 sprintf(Finley_ErrorMsg,"sample points of coefficient y don't match (%d,%d)",p.numQuad,elements->numElements);
83 }
84
85
86 /* check coefficient */
87
88 if (p.numEqu==1 && p.numComp==1) {
89 if (!isEmpty(d)) {
90 if (! isDataPointShapeEqual(d,0,dimensions)) {
91 Finley_ErrorCode=TYPE_ERROR;
92 sprintf(Finley_ErrorMsg,"coefficient d, rank 0 expected.");
93 }
94 }
95 if (!isEmpty(y)) {
96 if (! isDataPointShapeEqual(y,0,dimensions)) {
97 Finley_ErrorCode=TYPE_ERROR;
98 sprintf(Finley_ErrorMsg,"coefficient y, rank 0 expected.");
99 }
100 }
101 } else {
102 if (!isEmpty(d)) {
103 dimensions[0]=p.numEqu;
104 dimensions[1]=p.numComp;
105 if (! isDataPointShapeEqual(d,2,dimensions)) {
106 Finley_ErrorCode=TYPE_ERROR;
107 sprintf(Finley_ErrorMsg,"coefficient D, expected shape (%d,%d)",dimensions[0],dimensions[1]);
108 }
109 }
110 if (!isEmpty(y)) {
111 dimensions[0]=p.numEqu;
112 if (! isDataPointShapeEqual(y,1,dimensions)) {
113 Finley_ErrorCode=TYPE_ERROR;
114 sprintf(Finley_ErrorMsg,"coefficient D, expected shape (%d,)",dimensions[0]);
115 }
116 }
117 }
118
119 if (Finley_ErrorCode==NO_ERROR) {
120 time0=Finley_timer();
121 #pragma omp parallel private(index_row,index_col,EM_S,EM_F,V,dVdv,dSdV,Area,color)
122 {
123 EM_S=EM_F=V=dVdv=dSdV=Area=NULL;
124 index_row=index_col=NULL;
125 /* allocate work arrays: */
126 EM_S=THREAD_MEMALLOC(p.NN_col*p.NN_row*p.numEqu*p.numComp,double);
127 EM_F=THREAD_MEMALLOC(p.NN_row*p.numEqu,double);
128 V=THREAD_MEMALLOC(p.NN*p.numDim,double);
129 dVdv=THREAD_MEMALLOC(p.numDim*p.numElementDim*p.numQuad,double);
130 Area=THREAD_MEMALLOC(p.numQuad,double);
131 index_row=THREAD_MEMALLOC(p.NN_row,maybelong);
132 index_col=THREAD_MEMALLOC(p.NN_col,maybelong);
133 if (! ( Finley_checkPtr(EM_S) || Finley_checkPtr(EM_F) || Finley_checkPtr(index_row) || Finley_checkPtr(index_col) ||
134 Finley_checkPtr(V) || Finley_checkPtr(dVdv) || Finley_checkPtr(Area))) {
135 /* open loop over all colors: */
136 for (color=0;color<elements->numColors;color++) {
137 /* open loop over all elements: */
138 #pragma omp for private(e,q) schedule(static)
139 for(e=0;e<elements->numElements;e++){
140 if (elements->Color[e]==color) {
141 for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
142 /* gather V-coordinates of nodes into V: */
143 Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
144 /* handel the case where NS!=p.NN (contact elements) */
145 Finley_Assemble_handelShapeMissMatch_Mean_in(p.numDim,p.NN,1,V,p.NS,1);
146 /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
147 Finley_Util_SmallMatMult(p.numDim,p.numElementDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);
148 /* */
149 Finley_LengthOfNormalVector(p.numQuad,p.numDim,p.numElementDim,dVdv,Area);
150 /* scale area: */
151 for (q=0;q<p.numQuad;q++) Area[q]=ABS(Area[q]*p.referenceElement->QuadWeights[q]);
152 /* integration for the stiffness matrix: */
153 /* in order to optimze the number of operations the case of constants coefficience needs a bit more work */
154 /* to make use of some symmetry. */
155 if (S!=NULL) {
156 for (q=0;q<p.NN_col*p.NN_row*p.numEqu*p.numComp;q++) EM_S[q]=0;
157 if (p.numEqu==1 && p.numComp==1) {
158 Finley_Assemble_PDEMatrix_Single2(p.NS_row,p.numElementDim,p.numQuad,
159 p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_S,
160 NULL,TRUE,
161 NULL,TRUE,
162 NULL,TRUE,
163 getSampleData(d,e),isExpanded(d));
164 } else {
165 Finley_Assemble_PDEMatrix_System2(p.NS_row,p.numElementDim,p.numQuad,
166 p.numEqu,p.numComp,p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_S,
167 NULL,TRUE,
168 NULL,TRUE,
169 NULL,TRUE,
170 getSampleData(d,e),isExpanded(d));
171 }
172 /* handel the case of p.NS<NN */
173 handelShapeMissMatchForEM(p.numEqu*p.numComp,p.NN_row,p.NN_col,EM_S,p.NS_row,p.NS_col);
174 /*
175 {int g;
176 for(g=0;g<p.NN*p.NN*p.numEqu*p.numComp;g++) printf("%f ",EM_S[g]);
177 printf("\n");
178 }
179 */
180 /* add */
181 for (q=0;q<p.NN_col;q++) index_col[q]=p.label_col[elements->Nodes[INDEX2(p.col_node[q],e,p.NN)]];
182 Finley_SystemMatrix_add(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);
183 }
184 if (!isEmpty(F)) {
185 for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;
186 if (p.numEqu==1) {
187 Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numElementDim,p.numQuad,
188 p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_F,
189 NULL,TRUE,
190 getSampleData(y,e),isExpanded(y));
191 } else {
192 Finley_Assemble_RHSMatrix_System(p.NS_row,p.numElementDim,p.numQuad,
193 p.numEqu,p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_F,
194 NULL,TRUE,
195 getSampleData(y,e),isExpanded(y));
196 }
197 /* handel the case of NS<NN */
198 handelShapeMissMatchForEM(p.numEqu,p.NN_row,1,EM_F,p.NS_row,1);
199 /*
200 {int g;
201 for(g=0;g<p.NN*p.numEqu;g++) printf("%f ",EM_F[g]);
202 printf("\n");
203 }
204 */
205 /* add */
206 Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
207 }
208 }
209 } /* end of element loop: */
210 } /* end of color loop: */
211 }
212 THREAD_MEMFREE(EM_S);
213 THREAD_MEMFREE(EM_F);
214 THREAD_MEMFREE(V);
215 THREAD_MEMFREE(dVdv);
216 THREAD_MEMFREE(dSdV);
217 THREAD_MEMFREE(Area);
218 THREAD_MEMFREE(index_row);
219 THREAD_MEMFREE(index_col);
220 }
221 printf("timing: assemblage natural BC: %.4e sec\n",Finley_timer()-time0);
222 }
223 }
224 /*
225 * $Log$
226 * Revision 1.4 2004/12/15 07:08:32 jgs
227 * *** empty log message ***
228 *
229 *
230 *
231 */

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26