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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 471 - (hide annotations)
Fri Jan 27 01:33:02 2006 UTC (14 years, 1 month ago) by jgs
File MIME type: text/plain
File size: 8151 byte(s)
reorganise finley src tree to remove inc dir and src/finley directory

1 jgs 150 /*
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 jgs 82
15     /**************************************************************/
16    
17     /* assembles the system of numEq PDEs into the stiffness matrix S: */
18    
19     /* -div(A*grad u)-div(B*u)+C*grad u + D*u */
20    
21     /* -(A_{k,i,m,j} u_m,j)_i-(B_{k,i,m} u_m)_i+C_{k,m,j} u_m,j-D_{k,m} u_m */
22    
23     /* u has numComp components. */
24    
25     /* Shape of the coefficients: */
26    
27     /* A = numEqu x numDim x numComp x numDim */
28     /* B = numDim x numEqu x numComp */
29     /* C = numEqu x numDim x numComp */
30     /* D = numEqu x numComp */
31    
32    
33     /**************************************************************/
34    
35 jgs 150 /* Author: gross@access.edu.au */
36     /* Version: $Id$ */
37 jgs 82
38     /**************************************************************/
39    
40    
41     #include "Assemble.h"
42    
43     /**************************************************************/
44    
45 jgs 123 void Finley_Assemble_PDEMatrix_System2(dim_t NS,dim_t numDim,dim_t numQuad,dim_t numEqu,dim_t numComp,
46 jgs 82 double* S,double* DSDX, double* Vol,
47     int NN, double* EM_S,
48 jgs 123 double* A, bool_t extendedA,
49     double* B, bool_t extendedB,
50     double* C, bool_t extendedC,
51     double* D, bool_t extendedD ) {
52     dim_t s,r,k,m,i,j,q;
53 jgs 82 double rtmp;
54    
55    
56     /**************************************************************/
57     /* process A: */
58     /**************************************************************/
59     if (NULL!=A) {
60     if (extendedA) {
61     for (s=0;s<NS;s++) {
62     for (r=0;r<NS;r++) {
63     for (k=0;k<numEqu;k++) {
64     for (m=0;m<numComp;m++) {
65     for (i=0;i<numDim;i++) {
66     for (j=0;j<numDim;j++) {
67     for (q=0;q<numQuad;q++) {
68     EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=Vol[q]*DSDX[INDEX3(s,i,q,NS,numDim)]*
69     A[INDEX5(k,i,m,j,q,numEqu,numDim,numComp,numDim)]*DSDX[INDEX3(r,j,q,NS,numDim)];
70     }
71     }
72     }
73     }
74     }
75     }
76     }
77     } else {
78     for (s=0;s<NS;s++) {
79     for (r=0;r<NS;r++) {
80     for (i=0;i<numDim;i++) {
81     for (j=0;j<numDim;j++) {
82     rtmp=0;
83     for (q=0;q<numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,i,q,NS,numDim)]*DSDX[INDEX3(r,j,q,NS,numDim)];
84     for (k=0;k<numEqu;k++) {
85     for (m=0;m<numComp;m++) {
86     EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=rtmp*A[INDEX4(k,i,m,j,numEqu,numDim,numComp)];
87     }
88     }
89     }
90     }
91     }
92     }
93     }
94     }
95     /**************************************************************/
96     /* process B: */
97     /**************************************************************/
98     if (NULL!=B) {
99     if (extendedB) {
100     for (s=0;s<NS;s++) {
101     for (r=0;r<NS;r++) {
102     for (k=0;k<numEqu;k++) {
103     for (m=0;m<numComp;m++) {
104     for (i=0;i<numDim;i++) {
105     for (q=0;q<numQuad;q++) {
106     EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=Vol[q]*DSDX[INDEX3(s,i,q,NS,numDim)]*
107     B[INDEX4(k,i,m,q,numEqu,numDim,numComp)]*S[INDEX2(r,q,NS)];
108     }
109     }
110     }
111     }
112     }
113     }
114     } else {
115     for (s=0;s<NS;s++) {
116     for (r=0;r<NS;r++) {
117     for (i=0;i<numDim;i++) {
118     rtmp=0;
119     for (q=0;q<numQuad;q++) rtmp+=Vol[q]*DSDX[INDEX3(s,i,q,NS,numDim)]*S[INDEX2(r,q,NS)];
120     for (k=0;k<numEqu;k++) {
121     for (m=0;m<numComp;m++) {
122     EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=rtmp*B[INDEX3(k,i,m,numEqu,numDim)];
123     }
124     }
125     }
126     }
127     }
128     }
129     }
130     /**************************************************************/
131     /* process C: */
132     /**************************************************************/
133     if (NULL!=C) {
134     if (extendedC) {
135     for (s=0;s<NS;s++) {
136     for (r=0;r<NS;r++) {
137     for (k=0;k<numEqu;k++) {
138     for (m=0;m<numComp;m++) {
139     for (j=0;j<numDim;j++) {
140     for (q=0;q<numQuad;q++) {
141     EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=Vol[q]*S[INDEX2(s,q,NS)]*
142     C[INDEX4(k,m,j,q,numEqu,numComp,numDim)]*DSDX[INDEX3(r,j,q,NS,numDim)];
143     }
144     }
145     }
146     }
147     }
148     }
149     } else {
150     for (s=0;s<NS;s++) {
151     for (r=0;r<NS;r++) {
152     for (j=0;j<numDim;j++) {
153     rtmp=0;
154     for (q=0;q<numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,NS)]*DSDX[INDEX3(r,j,q,NS,numDim)];
155     for (k=0;k<numEqu;k++) {
156     for (m=0;m<numComp;m++) {
157     EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=rtmp*C[INDEX3(k,m,j,numEqu,numComp)];
158     }
159     }
160     }
161     }
162     }
163     }
164     }
165     /************************************************************* */
166     /* process D */
167     /**************************************************************/
168     if (NULL!=D) {
169     if (extendedD) {
170     for (s=0;s<NS;s++) {
171     for (r=0;r<NS;r++) {
172     for (k=0;k<numEqu;k++) {
173     for (m=0;m<numComp;m++) {
174     for (q=0;q<numQuad;q++) {
175     EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=Vol[q]*S[INDEX2(s,q,NS)]*
176     D[INDEX3(k,m,q,numEqu,numComp)]*S[INDEX2(r,q,NS)];
177     }
178     }
179     }
180     }
181     }
182     } else {
183     for (s=0;s<NS;s++) {
184     for (r=0;r<NS;r++) {
185     rtmp=0;
186     for (q=0;q<numQuad;q++) rtmp+=Vol[q]*S[INDEX2(s,q,NS)]*S[INDEX2(r,q,NS)];
187     for (k=0;k<numEqu;k++) {
188     for (m=0;m<numComp;m++) {
189     EM_S[INDEX4(k,m,s,r,numEqu,numComp,NN)]+=rtmp*D[INDEX2(k,m,numEqu)];
190     }
191     }
192     }
193     }
194     }
195     }
196     }
197    
198     /*
199     * $Log$
200 jgs 150 * Revision 1.3 2005/09/15 03:44:21 jgs
201     * Merge of development branch dev-02 back to main trunk on 2005-09-15
202     *
203     * Revision 1.2.2.1 2005/09/07 06:26:17 gross
204     * the solver from finley are put into the standalone package paso now
205     *
206 jgs 123 * Revision 1.2 2005/07/08 04:07:46 jgs
207     * Merge of development branch back to main trunk on 2005-07-08
208 jgs 82 *
209 jgs 123 * Revision 1.1.1.1.2.1 2005/06/29 02:34:47 gross
210     * some changes towards 64 integers in finley
211     *
212     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
213     * initial import of project esys2
214     *
215 jgs 82 * Revision 1.2 2004/07/30 04:37:06 gross
216     * escript and finley are linking now and RecMeshTest.py has been passed
217     *
218     * Revision 1.1 2004/07/02 04:21:13 gross
219     * Finley C code has been included
220     *
221     *
222     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26