/[escript]/trunk/esys2/finley/src/finleyC/Assemble_integrate.c
ViewVC logotype

Contents of /trunk/esys2/finley/src/finleyC/Assemble_integrate.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 100 - (show annotations)
Wed Dec 15 03:48:48 2004 UTC (14 years, 11 months ago) by jgs
File MIME type: text/plain
File size: 5482 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* assemblage routines: integrates data on quadrature points */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia, 2003,2004 */
10 /* author: gross@access.edu.au */
11 /* Version: $Id$ */
12
13 /**************************************************************/
14
15 #include "escript/Data/DataC.h"
16 #include "Finley.h"
17 #include "Assemble.h"
18 #include "NodeFile.h"
19 #include "ElementFile.h"
20 #include "Util.h"
21 #ifdef _OPENMP
22 #include <omp.h>
23 #endif
24
25 /**************************************************************/
26
27 void Finley_Assemble_integrate(Finley_NodeFile* nodes, Finley_ElementFile* elements,escriptDataC* data,double* out) {
28 int q,e,datacase,node_offset,i;
29 double *local_X, *dVdv, *Vol,*out_local,*data_array,rtmp;
30 if (nodes==NULL || elements==NULL) return;
31 int NN=elements->ReferenceElement->Type->numNodes;
32 int NS=elements->ReferenceElement->Type->numShapes;
33 int data_type=getFunctionSpaceType(data);
34 int numComps=getDataPointSize(data);
35 int numQuad=elements->ReferenceElement->numQuadNodes;
36 int numDim=nodes->numDim;
37 int numDim_local=elements->ReferenceElement->Type->numDim;
38
39 /* set some parameter */
40
41 if (data_type==FINLEY_ELEMENTS) {
42 datacase=0;
43 node_offset=0;
44 } else if (data_type==FINLEY_FACE_ELEMENTS) {
45 datacase=1;
46 node_offset=0;
47 } else if (data_type==FINLEY_CONTACT_ELEMENTS_1) {
48 datacase=1;
49 node_offset=0;
50 } else if (data_type==FINLEY_CONTACT_ELEMENTS_2) {
51 datacase=1;
52 node_offset=NN-NS;
53 } else {
54 Finley_ErrorCode=TYPE_ERROR;
55 sprintf(Finley_ErrorMsg,"integration of data is not possible.");
56 }
57
58 /* check the shape of the data */
59
60 if (! numSamplesEqual(data,numQuad,elements->numElements)) {
61 Finley_ErrorCode=TYPE_ERROR;
62 sprintf(Finley_ErrorMsg,"illegal number of samples of integrant kernel Data object");
63 }
64
65 /* now we can start */
66
67 if (Finley_ErrorCode==NO_ERROR) {
68 for (q=0;q<numComps;q++) out[q]=0;
69 #pragma omp parallel private(local_X,dVdv,Vol,out_local,i)
70 {
71
72 /* allocation of work arrays */
73
74 local_X=dVdv=Vol=out_local=NULL;
75 out_local=(double*) THREAD_MEMALLOC(numComps*sizeof(double));
76 dVdv=(double*) THREAD_MEMALLOC(numQuad*numDim_local*numDim*sizeof(double));
77 Vol=(double*) THREAD_MEMALLOC(numQuad*sizeof(double));
78 local_X=(double*) THREAD_MEMALLOC(NS*numDim*sizeof(double));
79
80 if (! (Finley_checkPtr(Vol) || Finley_checkPtr(local_X) || Finley_checkPtr(dVdv) || Finley_checkPtr(out_local))) {
81
82 /* initialize local result */
83
84 for (i=0;i<numComps;i++) out_local[i]=0;
85
86 /* open the element loop */
87
88 #pragma omp for private(e,q,data_array,rtmp) schedule(static)
89 for(e=0;e<elements->numElements;e++) {
90 /* gather local coordinates of nodes into local_X: */
91 Finley_Util_Gather_double(NS,&(elements->Nodes[INDEX2(node_offset,e,NN)]),numDim,nodes->Coordinates,local_X);
92 /* calculate dVdv(i,j,q)=V(i,k)*DSDv(k,j,q) */
93 Finley_Util_SmallMatMult(numDim,numDim_local*numQuad,dVdv,NS,local_X,elements->ReferenceElement->dSdv);
94 /* calculate volume /area elements */
95 switch (datacase) {
96 case 0:
97 Finley_Util_DetOfSmallMat(numQuad,numDim,dVdv,Vol);
98 break;
99 case 1:
100 Finley_LengthOfNormalVector(numQuad,numDim,numDim_local,dVdv,Vol);
101 break;
102 }
103 /* out=out+ integrals */
104 data_array=getSampleData(data,e);
105 if (isExpanded(data)) {
106 for (q=0;q<numQuad;q++) {
107 rtmp=ABS(Vol[q])*elements->ReferenceElement->QuadWeights[q];
108 for (i=0;i<numComps;i++) out_local[i]+=data_array[INDEX2(i,q,numComps)]*rtmp;
109 }
110 } else {
111 rtmp=0.;
112 for (q=0;q<numQuad;q++) rtmp+=ABS(Vol[q])*elements->ReferenceElement->QuadWeights[q];
113 for (i=0;i<numComps;i++) out_local[i]+=data_array[i]*rtmp;
114 }
115 }
116
117 /* add local results to global result */
118 #pragma omp critical
119 for (i=0;i<numComps;i++) out[i]+=out_local[i];
120
121 }
122 THREAD_MEMFREE(local_X);
123 THREAD_MEMFREE(dVdv);
124 THREAD_MEMFREE(Vol);
125 THREAD_MEMFREE(out_local);
126 }
127 }
128 }
129
130 /*
131 * $Log$
132 * Revision 1.3 2004/12/15 03:48:45 jgs
133 * *** empty log message ***
134 *
135 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
136 * initial import of project esys2
137 *
138 * Revision 1.2 2004/07/21 05:00:54 gross
139 * name changes in DataC
140 *
141 * Revision 1.1 2004/07/02 04:21:13 gross
142 * Finley C code has been included
143 *
144 *
145 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26