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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 155 - (hide annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 11 months 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 jgs 150 /*
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 jgs 82
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 jgs 150 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 jgs 82 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Area=NULL;
52     double time0;
53     Assemble_Parameters p;
54 jgs 123 index_t *index_row=NULL,*index_col=NULL,color;
55     dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
56 jgs 150 Finley_resetError();
57 jgs 82
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 jgs 150 if (! Finley_noError()) return;
64 jgs 82
65     /* get a functionspace */
66 jgs 123 type_t funcspace=UNKNOWN;
67 jgs 82 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 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient d");
75 jgs 82 }
76     if (! functionSpaceTypeEqual(funcspace,y) ) {
77 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient y");
78 jgs 82 }
79    
80     /* check if all function spaces are the same */
81    
82     if (! numSamplesEqual(d,p.numQuad,elements->numElements) ) {
83 jgs 150 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 jgs 82 }
86    
87     if (! numSamplesEqual(y,p.numQuad,elements->numElements) ) {
88 jgs 150 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 jgs 82 }
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 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: coefficient d, rank 0 expected.");
99 jgs 82 }
100     }
101     if (!isEmpty(y)) {
102     if (! isDataPointShapeEqual(y,0,dimensions)) {
103 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: coefficient y, rank 0 expected.");
104 jgs 82 }
105     }
106     } else {
107     if (!isEmpty(d)) {
108     dimensions[0]=p.numEqu;
109     dimensions[1]=p.numComp;
110     if (! isDataPointShapeEqual(d,2,dimensions)) {
111 jgs 150 sprintf(error_msg,"__FILE__: coefficient d, expected shape (%d,%d)",dimensions[0],dimensions[1]);
112     Finley_setError(TYPE_ERROR,error_msg);
113 jgs 82 }
114     }
115     if (!isEmpty(y)) {
116     dimensions[0]=p.numEqu;
117     if (! isDataPointShapeEqual(y,1,dimensions)) {
118 jgs 150 sprintf(error_msg,"__FILE__: coefficient y, expected shape (%d,)",dimensions[0]);
119     Finley_setError(TYPE_ERROR,error_msg);
120 jgs 82 }
121     }
122     }
123    
124 jgs 150 if (Finley_noError()) {
125 jgs 82 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 jgs 102 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 jgs 123 index_row=THREAD_MEMALLOC(p.NN_row,index_t);
137     index_col=THREAD_MEMALLOC(p.NN_col,index_t);
138 jgs 82 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 jgs 123 for (color=elements->minColor;color<=elements->maxColor;color++) {
142 jgs 82 /* 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 jgs 150 Finley_Assemble_addToSystemMatrix(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);
188 jgs 82 }
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 jgs 149 #ifdef Finley_TRACE
227 jgs 82 printf("timing: assemblage natural BC: %.4e sec\n",Finley_timer()-time0);
228 jgs 149 #endif
229 jgs 82 }
230     }
231     /*
232     * $Log$
233 jgs 150 * 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 jgs 149 * 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 jgs 147 * 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 jgs 150 * 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 jgs 149 * 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 jgs 147 * 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 jgs 123 * 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 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
258 jgs 97 * *** empty log message ***
259 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
260     * some changes towards 64 integers in finley
261 jgs 82 *
262 jgs 123 * 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 jgs 97 *
265 jgs 82 *
266 jgs 123 *
267 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26