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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 616 - (hide annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 6 months ago) by elspeth
File MIME type: text/plain
File size: 7852 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 147
13     /**************************************************************/
14    
15     /* assembles a system of numEq natural boundary condition into the stiffness matrix S and right hand side F: */
16    
17     /* y */
18    
19     /* Shape of the coefficients: */
20    
21     /* y = numEqu */
22    
23     /* The coefficient y have to be defined on the integartion points on face elements or not present (=NULL). */
24    
25     /* F have to be initialized before the routine is called.F can be NULL. */
26    
27    
28     /**************************************************************/
29    
30 jgs 150 /* Author: gross@access.edu.au */
31     /* Version: $Id$ */
32 jgs 147
33     /**************************************************************/
34    
35     #include "Assemble.h"
36     #include "Util.h"
37     #ifdef _OPENMP
38     #include <omp.h>
39     #endif
40    
41     /**************************************************************/
42    
43     void Finley_Assemble_RobinCondition_RHS(Finley_NodeFile* nodes,Finley_ElementFile* elements,escriptDataC* F, escriptDataC* y, Finley_Assemble_handelShapeMissMatch handelShapeMissMatchForEM) {
44 jgs 150 char error_msg[LenErrorMsg_MAX];
45 jgs 147 double *EM_F=NULL,*V=NULL,*dVdv=NULL,*dSdV=NULL,*Area=NULL;
46     double time0;
47     Assemble_Parameters p;
48     index_t *index_row=NULL,color;
49     dim_t dimensions[ESCRIPT_MAX_DATA_RANK],e,q;
50     bool_t assemble;
51 jgs 150 Finley_resetError();
52 jgs 147
53     if (nodes==NULL || elements==NULL) return;
54     if (isEmpty(F) || isEmpty(y) ) return;
55    
56     /* set all parameters in p*/
57     Assemble_getAssembleParameters(nodes,elements,NULL,F,&p);
58 jgs 150 if (! Finley_noError()) return;
59 jgs 147
60     /* get a functionspace */
61     type_t funcspace=UNKNOWN;
62     updateFunctionSpaceType(funcspace,y);
63     if (funcspace==UNKNOWN) return; /* all data are empty */
64    
65     /* check if all function spaces are the same */
66    
67     if (! functionSpaceTypeEqual(funcspace,y) ) {
68 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient y");
69 jgs 147 }
70    
71     /* check if all function spaces are the same */
72    
73     if (! numSamplesEqual(y,p.numQuad,elements->numElements) ) {
74 jgs 150 sprintf(error_msg,"__FILE__:sample points of coefficient y don't match (%d,%d)",p.numQuad,elements->numElements);
75     Finley_setError(TYPE_ERROR,error_msg);
76 jgs 147 }
77    
78    
79     /* check coefficient */
80    
81     if (p.numEqu==1 && p.numComp==1) {
82     if (!isEmpty(y)) {
83     if (! isDataPointShapeEqual(y,0,dimensions)) {
84 jgs 150 Finley_setError(TYPE_ERROR,"__FILE__: coefficient y, rank 0 expected.");
85 jgs 147 }
86     }
87     } else {
88     if (!isEmpty(y)) {
89     dimensions[0]=p.numEqu;
90     if (! isDataPointShapeEqual(y,1,dimensions)) {
91 jgs 150 sprintf(error_msg,"__FILE__:coefficient y, expected shape (%d,)",dimensions[0]);
92     Finley_setError(TYPE_ERROR,error_msg);
93 jgs 147 }
94     }
95     }
96    
97 jgs 150 if (Finley_noError()) {
98 jgs 147 time0=Finley_timer();
99     #pragma omp parallel private(assemble,index_row,EM_F,V,dVdv,dSdV,Area,color)
100     {
101     EM_F=V=dVdv=dSdV=Area=NULL;
102     index_row=NULL;
103     /* allocate work arrays: */
104     EM_F=THREAD_MEMALLOC(p.NN_row*p.numEqu,double);
105     V=THREAD_MEMALLOC(p.NN*p.numDim,double);
106     dVdv=THREAD_MEMALLOC(p.numDim*p.numElementDim*p.numQuad,double);
107     Area=THREAD_MEMALLOC(p.numQuad,double);
108     index_row=THREAD_MEMALLOC(p.NN_row,index_t);
109     if (! ( Finley_checkPtr(EM_F) || Finley_checkPtr(index_row) || Finley_checkPtr(V) || Finley_checkPtr(dVdv) || Finley_checkPtr(Area))) {
110     /* open loop over all colors: */
111     for (color=elements->minColor;color<=elements->maxColor;color++) {
112     /* open loop over all elements: */
113     #pragma omp for private(e,q) schedule(dynamic)
114     for(e=0;e<elements->numElements;e++){
115     if (elements->Color[e]==color) {
116     if (isExpanded(y)) {
117     assemble=Finley_Util_anyNonZeroDouble(p.numEqu*p.numQuad,getSampleData(y,e));
118     } else {
119     assemble=Finley_Util_anyNonZeroDouble(p.numEqu,getSampleData(y,e));
120     }
121     if (assemble) {
122     for (q=0;q<p.NN_row;q++) index_row[q]=p.label_row[elements->Nodes[INDEX2(p.row_node[q],e,p.NN)]];
123     /* gather V-coordinates of nodes into V: */
124     Finley_Util_Gather_double(p.NN,&(elements->Nodes[INDEX2(0,e,p.NN)]),p.numDim,nodes->Coordinates,V);
125     /* handel the case where NS!=p.NN (contact elements) */
126     Finley_Assemble_handelShapeMissMatch_Mean_in(p.numDim,p.NN,1,V,p.NS,1);
127     /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
128     Finley_Util_SmallMatMult(p.numDim,p.numElementDim*p.numQuad,dVdv,p.NS,V,p.referenceElement->dSdv);
129     /* */
130     Finley_LengthOfNormalVector(p.numQuad,p.numDim,p.numElementDim,dVdv,Area);
131     /* scale area: */
132     for (q=0;q<p.numQuad;q++) Area[q]=ABS(Area[q]*p.referenceElement->QuadWeights[q]);
133     for (q=0;q<p.NN_row*p.numEqu;q++) EM_F[q]=0;
134     if (p.numEqu==1) {
135     Finley_Assemble_RHSMatrix_Single(p.NS_row,p.numElementDim,p.numQuad,
136     p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_F,
137     NULL,TRUE,
138     getSampleData(y,e),isExpanded(y));
139     } else {
140     Finley_Assemble_RHSMatrix_System(p.NS_row,p.numElementDim,p.numQuad,
141     p.numEqu,p.referenceElement_row->S,dSdV,Area,p.NN_row,EM_F,
142     NULL,TRUE,
143     getSampleData(y,e),isExpanded(y));
144     }
145     /* handel the case of NS<NN */
146     handelShapeMissMatchForEM(p.numEqu,p.NN_row,1,EM_F,p.NS_row,1);
147     /* add */
148     Finley_Util_AddScatter(p.NN_row,index_row,p.numEqu,EM_F,getSampleData(F,0));
149     }
150     }
151     } /* end of element loop: */
152     } /* end of color loop: */
153     }
154     THREAD_MEMFREE(EM_F);
155     THREAD_MEMFREE(V);
156     THREAD_MEMFREE(dVdv);
157     THREAD_MEMFREE(dSdV);
158     THREAD_MEMFREE(Area);
159     THREAD_MEMFREE(index_row);
160     }
161 jgs 149 #ifdef Finley_TRACE
162 jgs 147 printf("timing: assemblage natural BC right hand side: %.4e sec\n",Finley_timer()-time0);
163 jgs 149 #endif
164 jgs 147 }
165     }
166     /*
167     * $Log$
168 jgs 150 * Revision 1.4 2005/09/15 03:44:21 jgs
169     * Merge of development branch dev-02 back to main trunk on 2005-09-15
170     *
171 jgs 149 * Revision 1.3 2005/09/01 03:31:35 jgs
172     * Merge of development branch dev-02 back to main trunk on 2005-09-01
173     *
174 jgs 147 * Revision 1.2 2005/08/12 01:45:43 jgs
175     * erge of development branch dev-02 back to main trunk on 2005-08-12
176     *
177 jgs 150 * Revision 1.1.2.3 2005/09/07 06:26:17 gross
178     * the solver from finley are put into the standalone package paso now
179     *
180 jgs 149 * Revision 1.1.2.2 2005/08/24 02:02:18 gross
181     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
182     *
183 jgs 147 * Revision 1.1.2.1 2005/08/04 22:41:11 gross
184     * some extra routines for finley that might speed-up RHS assembling in some cases (not actived right now)
185     *
186     *
187     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26