/[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 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (14 years, 11 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Assemble_RobinCondition.c
File MIME type: text/plain
File size: 10721 byte(s)
Initial revision

1 jgs 82 /* $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=(double*) THREAD_MEMALLOC(p.NN_col*p.NN_row*p.numEqu*p.numComp*sizeof(double));
127     EM_F=(double*) THREAD_MEMALLOC(p.NN_row*p.numEqu*sizeof(double));
128     V=(double*) THREAD_MEMALLOC(p.NN*p.numDim*sizeof(double));
129     dVdv=(double*) THREAD_MEMALLOC(p.numDim*p.numElementDim*p.numQuad*sizeof(double));
130     Area=(double*) THREAD_MEMALLOC(p.numQuad*sizeof(double));
131     index_row=(maybelong*) THREAD_MEMALLOC(p.NN_row*sizeof(maybelong));
132     index_col=(maybelong*) THREAD_MEMALLOC(p.NN_col*sizeof(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.1 2004/10/26 06:53:57 jgs
227     * Initial revision
228     *
229     * Revision 1.3 2004/07/30 04:37:06 gross
230     * escript and finley are linking now and RecMeshTest.py has been passed
231     *
232     * Revision 1.2 2004/07/21 05:00:54 gross
233     * name changes in DataC
234     *
235     * Revision 1.1 2004/07/02 04:21:13 gross
236     * Finley C code has been included
237     *
238     *
239     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26