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

1 jgs 150 /*
2 elspeth 616 ************************************************************
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 jgs 150 */
12 jgs 82
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 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) {
48     char error_msg[LenErrorMsg_MAX];
49 jgs 82 double *EM_S=NULL,*EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Area=NULL;
50     double time0;
51     Assemble_Parameters p;
52 jgs 123 index_t *index_row=NULL,*index_col=NULL,color;
53     dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
54 jgs 150 Finley_resetError();
55 jgs 82
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 jgs 150 if (! Finley_noError()) return;
62 jgs 82
63     /* get a functionspace */
64 jgs 123 type_t funcspace=UNKNOWN;
65 jgs 82 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 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient d");
73 jgs 82 }
74     if (! functionSpaceTypeEqual(funcspace,y) ) {
75 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient y");
76 jgs 82 }
77    
78     /* check if all function spaces are the same */
79    
80     if (! numSamplesEqual(d,p.numQuad,elements->numElements) ) {
81 jgs 150 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 jgs 82 }
84    
85     if (! numSamplesEqual(y,p.numQuad,elements->numElements) ) {
86 jgs 150 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 jgs 82 }
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 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: coefficient d, rank 0 expected.");
97 jgs 82 }
98     }
99     if (!isEmpty(y)) {
100     if (! isDataPointShapeEqual(y,0,dimensions)) {
101 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: coefficient y, rank 0 expected.");
102 jgs 82 }
103     }
104     } else {
105     if (!isEmpty(d)) {
106     dimensions[0]=p.numEqu;
107     dimensions[1]=p.numComp;
108     if (! isDataPointShapeEqual(d,2,dimensions)) {
109 jgs 150 sprintf(error_msg,"__FILE__: coefficient d, expected shape (%d,%d)",dimensions[0],dimensions[1]);
110     Finley_setError(TYPE_ERROR,error_msg);
111 jgs 82 }
112     }
113     if (!isEmpty(y)) {
114     dimensions[0]=p.numEqu;
115     if (! isDataPointShapeEqual(y,1,dimensions)) {
116 jgs 150 sprintf(error_msg,"__FILE__: coefficient y, expected shape (%d,)",dimensions[0]);
117     Finley_setError(TYPE_ERROR,error_msg);
118 jgs 82 }
119     }
120     }
121    
122 jgs 150 if (Finley_noError()) {
123 jgs 82 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 jgs 102 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 jgs 123 index_row=THREAD_MEMALLOC(p.NN_row,index_t);
135     index_col=THREAD_MEMALLOC(p.NN_col,index_t);
136 jgs 82 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 jgs 123 for (color=elements->minColor;color<=elements->maxColor;color++) {
140 jgs 82 /* 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 jgs 150 Finley_Assemble_addToSystemMatrix(S,p.NN_row,index_row,p.numEqu,p.NN_col,index_col,p.numComp,EM_S);
186 jgs 82 }
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 jgs 149 #ifdef Finley_TRACE
225 jgs 82 printf("timing: assemblage natural BC: %.4e sec\n",Finley_timer()-time0);
226 jgs 149 #endif
227 jgs 82 }
228     }
229     /*
230     * $Log$
231 jgs 150 * 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 jgs 149 * 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 jgs 147 * 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 jgs 150 * 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 jgs 149 * 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 jgs 147 * 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 jgs 123 * 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 jgs 102 * Revision 1.4 2004/12/15 07:08:32 jgs
256 jgs 97 * *** empty log message ***
257 jgs 123 * Revision 1.1.1.1.2.2 2005/06/29 02:34:47 gross
258     * some changes towards 64 integers in finley
259 jgs 82 *
260 jgs 123 * 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 jgs 97 *
263 jgs 82 *
264 jgs 123 *
265 jgs 82 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26