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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (14 years ago) by jgs
Original Path: trunk/finley/src/finleyC/Assemble_RobinCondition.c
File MIME type: text/plain
File size: 12314 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 /*
2 ******************************************************************************
3 * *
4 * COPYRIGHT ACcESS 2003,2004,2005 - All Rights Reserved *
5 * *
6 * This software is the property of ACcESS. No part of this code *
7 * may be copied in any form or by any means without the expressed written *
8 * consent of ACcESS. Copying, use or modification of this software *
9 * by any unauthorised person is illegal unless that person has a software *
10 * license agreement with ACcESS. *
11 * *
12 ******************************************************************************
13 */
14
15 /**************************************************************/
16
17 /* assembles a system of numEq natural boundary condition into the stiffness matrix S and right hand side F: */
18
19 /* d*u = y */
20
21 /* u has numComp components. */
22
23 /* Shape of the coefficients: */
24
25 /* d = numEqu x numComp */
26 /* y = numEqu */
27
28 /* The coefficients d and y have to be defined on the integartion points on face elements or not present (=NULL). */
29
30 /* S and F have to be initialized before the routine is called. S or F can be NULL. In this case the left or */
31 /* the right hand side of the natural boundary conditions is not processed. */
32
33
34 /**************************************************************/
35
36 /* author: gross@access.edu.au */
37 /* Version: $Id$ */
38
39 /**************************************************************/
40
41 #include "Assemble.h"
42 #include "Util.h"
43 #ifdef _OPENMP
44 #include <omp.h>
45 #endif
46
47 /**************************************************************/
48
49 void Finley_Assemble_RobinCondition(Finley_NodeFile* nodes,Finley_ElementFile* elements,Paso_SystemMatrix* S, escriptDataC* F, escriptDataC* d, escriptDataC* y, Finley_Assemble_handelShapeMissMatch handelShapeMissMatchForEM) {
50 char error_msg[LenErrorMsg_MAX];
51 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Area=NULL;
52 double time0;
53 Assemble_Parameters p;
54 index_t *index_row=NULL,*index_col=NULL,color;
55 dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
56 Finley_resetError();
57
58 if (nodes==NULL || elements==NULL) return;
59 if (S==NULL && isEmpty(F)) return;
60
61 /* set all parameters in p*/
62 Assemble_getAssembleParameters(nodes,elements,S,F,&p);
63 if (! Finley_noError()) return;
64
65 /* get a functionspace */
66 type_t funcspace=UNKNOWN;
67 updateFunctionSpaceType(funcspace,d);
68 updateFunctionSpaceType(funcspace,y);
69 if (funcspace==UNKNOWN) return; /* all data are empty */
70
71 /* check if all function spaces are the same */
72
73 if (! functionSpaceTypeEqual(funcspace,d) ) {
74 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient d");
75 }
76 if (! functionSpaceTypeEqual(funcspace,y) ) {
77 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient y");
78 }
79
80 /* check if all function spaces are the same */
81
82 if (! numSamplesEqual(d,p.numQuad,elements->numElements) ) {
83 sprintf(error_msg,"__FILE__: sample points of coefficient d don't match (%d,%d)",p.numQuad,elements->numElements);
84 Finley_setError(TYPE_ERROR,error_msg);
85 }
86
87 if (! numSamplesEqual(y,p.numQuad,elements->numElements) ) {
88 sprintf(error_msg,"__FILE__: sample points of coefficient y don't match (%d,%d)",p.numQuad,elements->numElements);
89 Finley_setError(TYPE_ERROR,error_msg);
90 }
91
92
93 /* check coefficient */
94
95 if (p.numEqu==1 && p.numComp==1) {
96 if (!isEmpty(d)) {
97 if (! isDataPointShapeEqual(d,0,dimensions)) {
98 Finley_setError(TYPE_ERROR,"__FILE__: coefficient d, rank 0 expected.");
99 }
100 }
101 if (!isEmpty(y)) {
102 if (! isDataPointShapeEqual(y,0,dimensions)) {
103 Finley_setError(TYPE_ERROR,"__FILE__: coefficient y, rank 0 expected.");
104 }
105 }
106 } else {
107 if (!isEmpty(d)) {
108 dimensions[0]=p.numEqu;
109 dimensions[1]=p.numComp;
110 if (! isDataPointShapeEqual(d,2,dimensions)) {
111 sprintf(error_msg,"__FILE__: coefficient d, expected shape (%d,%d)",dimensions[0],dimensions[1]);
112 Finley_setError(TYPE_ERROR,error_msg);
113 }
114 }
115 if (!isEmpty(y)) {
116 dimensions[0]=p.numEqu;
117 if (! isDataPointShapeEqual(y,1,dimensions)) {
118 sprintf(error_msg,"__FILE__: coefficient y, expected shape (%d,)",dimensions[0]);
119 Finley_setError(TYPE_ERROR,error_msg);
120 }
121 }
122 }
123
124 if (Finley_noError()) {
125 time0=Finley_timer();
126 #pragma omp parallel private(index_row,index_col,EM_S,EM_F,V,dVdv,dSdV,Area,color)
127 {
128 EM_S=EM_F=V=dVdv=dSdV=Area=NULL;
129 index_row=index_col=NULL;
130 /* allocate work arrays: */
131 EM_S=THREAD_MEMALLOC(p.NN_col*p.NN_row*p.numEqu*p.numComp,double);
132 EM_F=THREAD_MEMALLOC(p.NN_row*p.numEqu,double);
133 V=THREAD_MEMALLOC(p.NN*p.numDim,double);
134 dVdv=THREAD_MEMALLOC(p.numDim*p.numElementDim*p.numQuad,double);
135 Area=THREAD_MEMALLOC(p.numQuad,double);
136 index_row=THREAD_MEMALLOC(p.NN_row,index_t);
137 index_col=THREAD_MEMALLOC(p.NN_col,index_t);
138 if (! ( Finley_checkPtr(EM_S) || Finley_checkPtr(EM_F) || Finley_checkPtr(index_row) || Finley_checkPtr(index_col) ||
139 Finley_checkPtr(V) || Finley_checkPtr(dVdv) || Finley_checkPtr(Area))) {
140 /* open loop over all colors: */
141 for (color=elements->minColor;color<=elements->maxColor;color++) {
142 /* open loop over all elements: */
143 #pragma omp for private(e,q) schedule(static)
144 for(e=0;e<elements->numElements;e++){
145 if (elements->Color[e]==color) {
146 for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
147 /* gather V-coordinates of nodes into V: */
148 Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
149 /* handel the case where NS!=p.NN (contact elements) */
150 Finley_Assemble_handelShapeMissMatch_Mean_in(p.numDim,p.NN,1,V,p.NS,1);
151 /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
152 Finley_Util_SmallMatMult(p.numDim,p.numElementDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);
153 /* */
154 Finley_LengthOfNormalVector(p.numQuad,p.numDim,p.numElementDim,dVdv,Area);
155 /* scale area: */
156 for (q=0;q<p.numQuad;q++) Area[q]=ABS(Area[q]*p.referenceElement->QuadWeights[q]);
157 /* integration for the stiffness matrix: */
158 /* in order to optimze the number of operations the case of constants coefficience needs a bit more work */
159 /* to make use of some symmetry. */
160 if (S!=NULL) {
161 for (q=0;q<p.NN_col*p.NN_row*p.numEqu*p.numComp;q++) EM_S[q]=0;
162 if (p.numEqu==1 && p.numComp==1) {
163 Finley_Assemble_PDEMatrix_Single2(p.NS_row,p.numElementDim,p.numQuad,
164 p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_S,
165 NULL,TRUE,
166 NULL,TRUE,
167 NULL,TRUE,
168 getSampleData(d,e),isExpanded(d));
169 } else {
170 Finley_Assemble_PDEMatrix_System2(p.NS_row,p.numElementDim,p.numQuad,
171 p.numEqu,p.numComp,p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_S,
172 NULL,TRUE,
173 NULL,TRUE,
174 NULL,TRUE,
175 getSampleData(d,e),isExpanded(d));
176 }
177 /* handel the case of p.NS<NN */
178 handelShapeMissMatchForEM(p.numEqu*p.numComp,p.NN_row,p.NN_col,EM_S,p.NS_row,p.NS_col);
179 /*
180 {int g;
181 for(g=0;g<p.NN*p.NN*p.numEqu*p.numComp;g++) printf("%f ",EM_S[g]);
182 printf("\n");
183 }
184 */
185 /* add */
186 for (q=0;q<p.NN_col;q++) index_col[q]=p.label_col[elements->Nodes[INDEX2(p.col_node[q],e,p.NN)]];
187 Finley_Assemble_addToSystemMatrix(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);
188 }
189 if (!isEmpty(F)) {
190 for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;
191 if (p.numEqu==1) {
192 Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numElementDim,p.numQuad,
193 p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_F,
194 NULL,TRUE,
195 getSampleData(y,e),isExpanded(y));
196 } else {
197 Finley_Assemble_RHSMatrix_System(p.NS_row,p.numElementDim,p.numQuad,
198 p.numEqu,p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_F,
199 NULL,TRUE,
200 getSampleData(y,e),isExpanded(y));
201 }
202 /* handel the case of NS<NN */
203 handelShapeMissMatchForEM(p.numEqu,p.NN_row,1,EM_F,p.NS_row,1);
204 /*
205 {int g;
206 for(g=0;g<p.NN*p.numEqu;g++) printf("%f ",EM_F[g]);
207 printf("\n");
208 }
209 */
210 /* add */
211 Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
212 }
213 }
214 } /* end of element loop: */
215 } /* end of color loop: */
216 }
217 THREAD_MEMFREE(EM_S);
218 THREAD_MEMFREE(EM_F);
219 THREAD_MEMFREE(V);
220 THREAD_MEMFREE(dVdv);
221 THREAD_MEMFREE(dSdV);
222 THREAD_MEMFREE(Area);
223 THREAD_MEMFREE(index_row);
224 THREAD_MEMFREE(index_col);
225 }
226 #ifdef Finley_TRACE
227 printf("timing: assemblage natural BC: %.4e sec\n",Finley_timer()-time0);
228 #endif
229 }
230 }
231 /*
232 * $Log$
233 * Revision 1.8 2005/09/15 03:44:21 jgs
234 * Merge of development branch dev-02 back to main trunk on 2005-09-15
235 *
236 * Revision 1.7 2005/09/01 03:31:35 jgs
237 * Merge of development branch dev-02 back to main trunk on 2005-09-01
238 *
239 * Revision 1.6 2005/08/12 01:45:43 jgs
240 * erge of development branch dev-02 back to main trunk on 2005-08-12
241 *
242 * Revision 1.5.2.4 2005/09/07 06:26:17 gross
243 * the solver from finley are put into the standalone package paso now
244 *
245 * Revision 1.5.2.3 2005/08/24 02:02:18 gross
246 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
247 *
248 * Revision 1.5.2.2 2005/08/04 22:41:11 gross
249 * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now)
250 *
251 * Revision 1.5.2.1 2005/08/03 08:55:46 gross
252 * typo fixed
253 *
254 * Revision 1.5 2005/07/08 04:07:46 jgs
255 * Merge of development branch back to main trunk on 2005-07-08
256 *
257 * Revision 1.4 2004/12/15 07:08:32 jgs
258 * *** empty log message ***
259 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
260 * some changes towards 64 integers in finley
261 *
262 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
263 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
264 *
265 *
266 *
267 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26