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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 97 - (hide annotations)
Tue Dec 14 05:39:33 2004 UTC (14 years, 11 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Assemble_integrate.c
File MIME type: text/plain
File size: 5567 byte(s)
*** empty log message ***

1 jgs 82 /* $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 jgs 97 out_local=THREAD_MEMALLOC(numComps,double);
76     dVdv=THREAD_MEMALLOC(numQuad*numDim_local*numDim,double);
77     Vol=THREAD_MEMALLOC(numQuad,double);
78     local_X=THREAD_MEMALLOC(NS*numDim,double);
79 jgs 82
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 jgs 97 * Revision 1.2 2004/12/14 05:39:30 jgs
133     * *** empty log message ***
134 jgs 82 *
135 jgs 97 * Revision 1.1.1.1.2.1 2004/11/24 01:37:12 gross
136     * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
137     *
138     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
139     * initial import of project esys2
140     *
141 jgs 82 * Revision 1.2 2004/07/21 05:00:54 gross
142     * name changes in DataC
143     *
144     * Revision 1.1 2004/07/02 04:21:13 gross
145     * Finley C code has been included
146     *
147     *
148     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26