/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Assemble_RobinCondition_RHS.c
File MIME type: text/plain
File size: 8174 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26