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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 201 - (show annotations)
Wed Nov 23 04:10:21 2005 UTC (14 years, 4 months ago) by jgs
File MIME type: text/plain
File size: 8174 byte(s)
copy finleyC and CPPAdapter to finley and finley/CPPAdapter to
facilitate scons builds

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