/[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 616 - (show annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 5 months ago) by elspeth
File MIME type: text/plain
File size: 11992 byte(s)
Copyright added to more source files.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26