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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 616 - (show 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 /*
2 ************************************************************
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 */
12
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 /* Author: gross@access.edu.au */
31 /* Version: $Id$ */
32
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 char error_msg[LenErrorMsg_MAX];
45 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 Finley_resetError();
52
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 if (! Finley_noError()) return;
59
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 Finley_setError(TYPE_ERROR,"__FILE__: unexpected function space type for coefficient y");
69 }
70
71 /* check if all function spaces are the same */
72
73 if (! numSamplesEqual(y,p.numQuad,elements->numElements) ) {
74 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 }
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 Finley_setError(TYPE_ERROR,"__FILE__: coefficient y, rank 0 expected.");
85 }
86 }
87 } else {
88 if (!isEmpty(y)) {
89 dimensions[0]=p.numEqu;
90 if (! isDataPointShapeEqual(y,1,dimensions)) {
91 sprintf(error_msg,"__FILE__:coefficient y, expected shape (%d,)",dimensions[0]);
92 Finley_setError(TYPE_ERROR,error_msg);
93 }
94 }
95 }
96
97 if (Finley_noError()) {
98 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 #ifdef Finley_TRACE
162 printf("timing: assemblage natural BC right hand side: %.4e sec\n",Finley_timer()-time0);
163 #endif
164 }
165 }
166 /*
167 * $Log$
168 * 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 * 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 * 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 * 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 * 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 * 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