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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (14 years, 6 months ago) by jgs
File MIME type: text/plain
File size: 5353 byte(s)
*** empty log message ***

1 /* $Id$ */
2
3 /**************************************************************/
4
5 /* Finley: Mesh */
6
7 /**************************************************************/
8
9 /* Copyrights by ACcESS Australia 2003, 2004 */
10 /* Author: gross@access.edu.au */
11
12 /**************************************************************/
13
14 #include "Finley.h"
15 #include "IndexList.h"
16 #include "System.h"
17 #include "Mesh.h"
18
19 /**************************************************************/
20
21 /* returns a reference to the matrix pattern */
22
23 Finley_SystemMatrixPattern* Finley_getPattern(Finley_Mesh *mesh,int reduce_row_order, int reduce_col_order) {
24 Finley_SystemMatrixPattern *out=NULL;
25 Finley_ErrorCode=NO_ERROR;
26 /* make sure that the requested pattern is available */
27 if (reduce_row_order) {
28 if (reduce_col_order) {
29 if (mesh->ReducedReducedPattern==NULL) mesh->ReducedReducedPattern=Finley_makePattern(mesh,reduce_row_order,reduce_col_order);
30 } else {
31 if (mesh->ReducedFullPattern==NULL) mesh->ReducedFullPattern=Finley_makePattern(mesh,reduce_row_order,reduce_col_order);
32 }
33 } else {
34 if (reduce_col_order) {
35 if (mesh->FullReducedPattern==NULL) mesh->FullReducedPattern=Finley_makePattern(mesh,reduce_row_order,reduce_col_order);
36 } else {
37 if (mesh->FullFullPattern==NULL) mesh->FullFullPattern=Finley_makePattern(mesh,reduce_row_order,reduce_col_order);
38 }
39 }
40 if (Finley_ErrorCode==NO_ERROR) {
41 if (reduce_row_order) {
42 if (reduce_col_order) {
43 out=Finley_SystemMatrixPattern_reference(mesh->ReducedReducedPattern);
44 } else {
45 out=Finley_SystemMatrixPattern_reference(mesh->ReducedFullPattern);
46 }
47 } else {
48 if (reduce_col_order) {
49 out=Finley_SystemMatrixPattern_reference(mesh->FullReducedPattern);
50 } else {
51 out=Finley_SystemMatrixPattern_reference(mesh->FullFullPattern);
52 }
53 }
54 }
55 return out;
56 }
57 Finley_SystemMatrixPattern* Finley_makePattern(Finley_Mesh *mesh,int reduce_row_order, int reduce_col_order) {
58 double time0;
59 maybelong i,n,s,itmp;
60 maybelong *rowLabel=NULL,*colLabel=NULL, *ptr=NULL, *index=NULL;
61 Finley_IndexList* index_list=NULL;
62 Finley_ErrorCode=NO_ERROR;
63
64 time0=Finley_timer();
65 if (reduce_col_order) {
66 n=mesh->Nodes->reducedNumDegreesOfFreedom;
67 colLabel=mesh->Nodes->reducedDegreeOfFreedom;
68 } else {
69 n=mesh->Nodes->numDegreesOfFreedom;
70 colLabel=mesh->Nodes->degreeOfFreedom;
71 }
72
73 if (reduce_row_order) {
74 n=mesh->Nodes->reducedNumDegreesOfFreedom;
75 rowLabel=mesh->Nodes->reducedDegreeOfFreedom;
76 } else {
77 n=mesh->Nodes->numDegreesOfFreedom;
78 rowLabel=mesh->Nodes->degreeOfFreedom;
79 }
80
81 index_list=TMPMEMALLOC(n,Finley_IndexList);
82 ptr=MEMALLOC(n+1,maybelong);
83 if (! (Finley_checkPtr(index_list) || Finley_checkPtr(ptr)) ) {
84 #pragma omp parallel private(i,s,itmp)
85 {
86 #pragma omp for schedule(static)
87 for(i=0;i<n;++i) {
88 index_list[i].extension=NULL;
89 index_list[i].n=0;
90 }
91 /* insert contributions from element matrices into colums index index_list: */
92 Finley_IndexList_insertElements(index_list,mesh->Elements,
93 reduce_row_order,rowLabel,reduce_col_order,colLabel);
94 Finley_IndexList_insertElements(index_list,mesh->FaceElements,
95 reduce_row_order,rowLabel,reduce_col_order,colLabel);
96 Finley_IndexList_insertElements(index_list,mesh->ContactElements,
97 reduce_row_order,rowLabel,reduce_col_order,colLabel);
98 Finley_IndexList_insertElements(index_list,mesh->Points,
99 reduce_row_order,rowLabel,reduce_col_order,colLabel);
100 /* get the number of connections per row */
101 #pragma omp for schedule(static)
102 for(i=0;i<n;++i) ptr[i]=Finley_IndexList_count(&index_list[i]);
103 /* accumalate ptr */
104 /* OMP */
105 #pragma omp master
106 {
107 s=0;
108 for(i=0;i<n;++i) {
109 itmp=ptr[i];
110 ptr[i]=s;
111 s+=itmp;
112 }
113 ptr[n]=s;
114 /* fill index */
115 index=MEMALLOC(ptr[n],maybelong);
116 }
117 #pragma omp barrier
118 if (! Finley_checkPtr(index)) {
119 #pragma omp for schedule(static)
120 for(i=0;i<n;++i) Finley_IndexList_toArray(&index_list[i],&index[ptr[i]]);
121 }
122 }
123 }
124 /* all done */
125 /* clean up */
126 if (index_list!=NULL) {
127 #pragma omp parallel for private(i)
128 for(i=0;i<n;++i) Finley_IndexList_free(index_list[i].extension);
129 }
130 TMPMEMFREE(index_list);
131 printf("timing: mesh to matrix pattern: %.4e sec\n",Finley_timer()-time0);
132 if (Finley_ErrorCode!=NO_ERROR) {
133 MEMFREE(ptr);
134 MEMFREE(index);
135 return NULL;
136 } else {
137 return Finley_SystemMatrixPattern_alloc(n,ptr,index);
138 }
139 }
140 /*
141 * $Log$
142 * Revision 1.2 2004/12/15 07:08:33 jgs
143 * *** empty log message ***
144 *
145 * Revision 1.1.2.3 2004/11/24 01:37:14 gross
146 * some changes dealing with the integer overflow in memory allocation. Finley solves 4M unknowns now
147 *
148 * Revision 1.1.2.1 2004/11/15 01:01:09 gross
149 * and anotther missing file
150 *
151 */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26