/[escript]/branches/arrayview_from_1695_trunk/finley/src/Assemble_PDE_System2_C.c
ViewVC logotype

Annotation of /branches/arrayview_from_1695_trunk/finley/src/Assemble_PDE_System2_C.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1696 - (hide annotations)
Mon Aug 11 03:33:40 2008 UTC (11 years ago) by jfenwick
File MIME type: text/plain
File size: 7818 byte(s)
Branching to perform experiments on removing DataArrayView

1 gross 798
2 ksteube 1312 /* $Id$ */
3    
4     /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16 gross 798 /**************************************************************/
17    
18     /* assembles the system of numEq PDEs into the stiffness matrix S right hand side F */
19    
20     /* D_{k,m} u_m and Y_k */
21    
22     /* u has p.numComp components in a 3D domain. The shape functions for test and solution must be identical */
23     /* and 2* row_NS == row_NN (contact elements) */
24    
25     /* Shape of the coefficients: */
26    
27     /* D = p.numEqu x p.numComp */
28     /* Y = p.numEqu */
29    
30    
31     /**************************************************************/
32    
33     /* Author: gross@access.edu.au */
34     /* Version: $Id:$ */
35    
36     /**************************************************************/
37    
38    
39     #include "Assemble.h"
40     #include "Util.h"
41 gross 853 #ifdef _OPENMP
42     #include <omp.h>
43     #endif
44 gross 798
45 gross 853
46 gross 798 /**************************************************************/
47    
48     void Finley_Assemble_PDE_System2_C(Assemble_Parameters p, Finley_ElementFile* elements,
49     Paso_SystemMatrix* Mat, escriptDataC* F, escriptDataC* D, escriptDataC* Y) {
50    
51     index_t color;
52     dim_t e;
53 gross 853 double *EM_S, *EM_F, *Vol, *D_p, *Y_p;
54     index_t *row_index;
55     register dim_t q, s,r,k,m;
56     register double rtmp, rtmp_D;
57     bool_t add_EM_F, add_EM_S;
58    
59 gross 798 bool_t extendedD=isExpanded(D);
60     bool_t extendedY=isExpanded(Y);
61     double *F_p=getSampleData(F,0);
62     double *S=p.row_jac->ReferenceElement->S;
63    
64    
65 gross 853 #pragma omp parallel private(color,EM_S, EM_F, Vol, D_p, Y_p,row_index,q, s,r,k,m,rtmp, rtmp_D,add_EM_F, add_EM_S)
66 gross 798 {
67 gross 853 EM_S=THREAD_MEMALLOC(p.row_NN*p.col_NN*p.numEqu*p.numComp,double);
68     EM_F=THREAD_MEMALLOC(p.row_NN*p.numEqu,double);
69     row_index=THREAD_MEMALLOC(p.row_NN,index_t);
70    
71     if (!Finley_checkPtr(EM_S) && !Finley_checkPtr(EM_F) && !Finley_checkPtr(row_index) ) {
72    
73     for (color=elements->minColor;color<=elements->maxColor;color++) {
74     /* open loop over all elements: */
75     #pragma omp for private(e) schedule(static)
76     for(e=0;e<elements->numElements;e++){
77     if (elements->Color[e]==color) {
78     Vol=&(p.row_jac->volume[INDEX2(0,e,p.numQuad)]);
79     add_EM_F=FALSE;
80     add_EM_S=FALSE;
81     /************************************************************* */
82     /* process D */
83     /**************************************************************/
84     D_p=getSampleData(D,e);
85     if (NULL!=D_p) {
86     add_EM_S=TRUE;
87     if (extendedD) {
88     for (s=0;s<p.row_NS;s++) {
89     for (r=0;r<p.col_NS;r++) {
90     for (k=0;k<p.numEqu;k++) {
91     for (m=0;m<p.numComp;m++) {
92     rtmp=0;
93     for (q=0;q<p.numQuad;q++) {
94     rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*D_p[INDEX3(k,m,q,p.numEqu,p.numComp)]*S[INDEX2(r,q,p.row_NS)];
95     }
96     EM_S[INDEX4(k,m,s ,r ,p.numEqu,p.numComp,p.row_NN)]= rtmp;
97     EM_S[INDEX4(k,m,s ,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]=-rtmp;
98     EM_S[INDEX4(k,m,s+p.row_NS,r ,p.numEqu,p.numComp,p.row_NN)]=-rtmp;
99     EM_S[INDEX4(k,m,s+p.row_NS,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]= rtmp;
100     }
101     }
102     }
103     }
104     } else {
105     for (s=0;s<p.row_NS;s++) {
106     for (r=0;r<p.col_NS;r++) {
107     rtmp=0;
108     for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*S[INDEX2(r,q,p.row_NS)];
109     for (k=0;k<p.numEqu;k++) {
110     for (m=0;m<p.numComp;m++) {
111     rtmp_D=rtmp*D_p[INDEX2(k,m,p.numEqu)];
112     EM_S[INDEX4(k,m,s ,r ,p.numEqu,p.numComp,p.row_NN)]= rtmp_D;
113     EM_S[INDEX4(k,m,s ,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]=-rtmp_D;
114     EM_S[INDEX4(k,m,s+p.row_NS,r ,p.numEqu,p.numComp,p.row_NN)]=-rtmp_D;
115     EM_S[INDEX4(k,m,s+p.row_NS,r+p.col_NS,p.numEqu,p.numComp,p.row_NN)]= rtmp_D;
116     }
117     }
118     }
119     }
120     }
121     }
122     /**************************************************************/
123     /* process Y: */
124     /**************************************************************/
125     Y_p=getSampleData(Y,e);
126     if (NULL!=Y_p) {
127     add_EM_F=TRUE;
128     if (extendedY) {
129     for (s=0;s<p.row_NS;s++) {
130     for (k=0;k<p.numEqu;k++) {
131 gross 798 rtmp=0;
132 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)]*Y_p[INDEX2(k,q,p.numEqu)];
133     EM_F[INDEX2(k,s ,p.numEqu)]=-rtmp;
134     EM_F[INDEX2(k,s+p.row_NS,p.numEqu)]= rtmp;
135     }
136 gross 798 }
137 gross 853 } else {
138     for (s=0;s<p.row_NS;s++) {
139 gross 798 rtmp=0;
140 gross 853 for (q=0;q<p.numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,p.row_NS)];
141 gross 798 for (k=0;k<p.numEqu;k++) {
142 gross 853 rtmp_D=rtmp*Y_p[k];
143     EM_F[INDEX2(k,s ,p.numEqu)]=-rtmp_D;
144     EM_F[INDEX2(k,s+p.row_NS,p.numEqu)]= rtmp_D;
145 gross 798 }
146     }
147     }
148 gross 853 }
149     /***********************************************************************************************/
150     /* add the element matrices onto the matrix and right hand side */
151     /***********************************************************************************************/
152     for (q=0;q<p.row_NN;q++) row_index[q]=p.row_DOF[elements->Nodes[INDEX2(p.row_node[q],e,p.row_NN)]];
153     if (add_EM_F) Finley_Util_AddScatter(p.row_NN,row_index,p.numEqu,EM_F,F_p, p.row_DOF_UpperBound);
154     if (add_EM_S) Finley_Assemble_addToSystemMatrix(Mat,p.row_NN,row_index,p.numEqu,p.col_NN,row_index,p.numComp,EM_S);
155    
156     } /* end color check */
157     } /* end element loop */
158     } /* end color loop */
159    
160     THREAD_MEMFREE(EM_S);
161     THREAD_MEMFREE(EM_F);
162     THREAD_MEMFREE(row_index);
163 gross 798
164 gross 853 } /* end of pointer check */
165 gross 798 } /* end parallel region */
166     }
167     /*
168     * $Log$
169     */

  ViewVC Help
Powered by ViewVC 1.1.26