SPEED
MAKE_SPX_CON.f90
Go to the documentation of this file.
1! Copyright (C) 2012 The SPEED FOUNDATION
2! Author: Ilario Mazzieri
3!
4! This file is part of SPEED.
5!
6! SPEED is free software; you can redistribute it and/or modify it
7! under the terms of the GNU Affero General Public License as
8! published by the Free Software Foundation, either version 3 of the
9! License, or (at your option) any later version.
10!
11! SPEED is distributed in the hope that it will be useful, but
12! WITHOUT ANY WARRANTY; without even the implied warranty of
13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14! Affero General Public License for more details.
15!
16! You should have received a copy of the GNU Affero General Public License
17! along with SPEED. If not, see <http://www.gnu.org/licenses/>.
18
30
31
32 subroutine make_spx_con(nelem,con_mac,nmat,tag_mat,sdeg,&
33 nnz_pntr,node_pntr,con_nnz,con_spx,nnode)
34
35 implicit none
36
37 integer*4 :: nelem,nmat,nnz_pntr,con_nnz,nnode
38 integer*4, dimension(nelem,9) :: con_mac
39 integer*4, dimension(nmat) :: tag_mat
40 integer*4, dimension(nmat) :: sdeg
41 integer*4, dimension(0:nnz_pntr) :: node_pntr
42 integer*4, dimension(0:con_nnz) :: con_spx
43
44 integer*4 :: nnode_mac,nn,imat,ie,je
45 integer*4 :: i,j,k,an,bn,cn,ia,ib,ic,ja,jb,jc,ka,kb,kc
46 integer*4 :: jsa,dljs,dmjs,l,m,is,js
47
48
49 do i = 0,con_nnz
50 con_spx(i) = 0
51 enddo
52
53 con_spx(0) = nelem +1
54 do ie = 1,nelem
55 do imat = 1,nmat
56 if (con_mac(ie,1).eq.tag_mat(imat)) then
57 nn = sdeg(imat) +1
58 con_spx(ie) = con_spx(ie -1) + nn*nn*nn +1
59 endif
60 enddo
61 enddo
62
63 nnode_mac = node_pntr(0) -1
64
65 nnode = nnode_mac
66
67 do ie = 1,nelem
68 do imat = 1,nmat
69 if (con_mac(ie,1).eq.tag_mat(imat)) then
70 nn = sdeg(imat) +1
71
72! First put material index
73 con_spx(con_spx(ie -1) +0) = con_mac(ie,1)
74
75! Then put vertices
76 con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +1) = &
77 con_mac(ie,2)
78 con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +nn) = &
79 con_mac(ie,3)
80 con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +nn) = &
81 con_mac(ie,4)
82 con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +1) = &
83 con_mac(ie,5)
84 con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +1) = &
85 con_mac(ie,6)
86 con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +nn) = &
87 con_mac(ie,7)
88 con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +nn) = &
89 con_mac(ie,8)
90 con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +1) = &
91 con_mac(ie,9)
92
93
94! Then construct edge connectivity
95
96! First edge down
97
98 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +1)
99 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +nn)
100
101 je = ie
102 do i = node_pntr(an -1),node_pntr(an) -1
103 do j = node_pntr(bn -1),node_pntr(bn) -1
104 if (node_pntr(i).eq.node_pntr(j)) then
105 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
106 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
107 je = node_pntr(i)
108 endif
109 endif
110 enddo
111 enddo
112
113 if (je.ne.ie) then
114 ia = 0
115 ja = 0
116 ka = 0
117 ib = 0
118 jb = 0
119 kb = 0
120
121 do k = 1,nn
122 do j = 1,nn
123 do i = 1,nn
124 js = nn*nn*(k -1) +nn*(j -1) +i
125
126 if (con_spx(con_spx(je -1) +js).eq.an) then
127 ia = i
128 ja = j
129 ka = k
130 endif
131
132 if (con_spx(con_spx(je -1) +js).eq.bn) then
133 ib = i
134 jb = j
135 kb = k
136 endif
137 enddo
138 enddo
139 enddo
140
141 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
142 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
143
144 do l = 1,nn
145 is = nn*nn*(1 -1) +nn*(1 -1) +l
146 js = jsa +dljs*(l -1)
147 con_spx(con_spx(ie -1) +is) = &
148 con_spx(con_spx(je -1) +js)
149 enddo
150
151 else
152 do i = 2,(nn -1)
153 nnode = nnode +1
154 is = nn*nn*(1 -1) +nn*(1 -1) +i
155 con_spx(con_spx(ie -1) + is) = nnode
156 enddo
157 endif
158
159
160! Second edge down
161
162 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +nn)
163 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +nn)
164
165 je = ie
166 do i = node_pntr(an -1),node_pntr(an) -1
167 do j = node_pntr(bn -1),node_pntr(bn) -1
168 if (node_pntr(i).eq.node_pntr(j)) then
169 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
170 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
171 je = node_pntr(i)
172 endif
173 endif
174 enddo
175 enddo
176
177 if (je.ne.ie) then
178 ia = 0
179 ja = 0
180 ka = 0
181 ib = 0
182 jb = 0
183 kb = 0
184
185 do k = 1,nn
186 do j = 1,nn
187 do i = 1,nn
188 js = nn*nn*(k -1) +nn*(j -1) +i
189
190 if (con_spx(con_spx(je -1) +js).eq.an) then
191 ia = i
192 ja = j
193 ka = k
194 endif
195
196 if (con_spx(con_spx(je -1) +js).eq.bn) then
197 ib = i
198 jb = j
199 kb = k
200 endif
201 enddo
202 enddo
203 enddo
204
205 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
206 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
207
208 do l = 1,nn
209 is = nn*nn*(1 -1) +nn*(l -1) +nn
210 js = jsa +dljs*(l -1)
211 con_spx(con_spx(ie -1) +is) = &
212 con_spx(con_spx(je -1) +js)
213 enddo
214
215 else
216 do j = 2,(nn -1)
217 nnode = nnode +1
218 is = nn*nn*(1 -1) +nn*(j -1) +nn
219 con_spx(con_spx(ie -1) + is) = nnode
220 enddo
221 endif
222
223
224! Third edge down
225
226 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +1)
227 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +nn)
228
229 je = ie
230 do i = node_pntr(an -1),node_pntr(an) -1
231 do j = node_pntr(bn -1),node_pntr(bn) -1
232 if (node_pntr(i).eq.node_pntr(j)) then
233 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
234 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
235 je = node_pntr(i)
236 endif
237 endif
238 enddo
239 enddo
240
241
242 if (je.ne.ie) then
243 ia = 0
244 ja = 0
245 ka = 0
246 ib = 0
247 jb = 0
248 kb = 0
249
250 do k = 1,nn
251 do j = 1,nn
252 do i = 1,nn
253 js = nn*nn*(k -1) +nn*(j -1) +i
254
255 if (con_spx(con_spx(je -1) +js).eq.an) then
256 ia = i
257 ja = j
258 ka = k
259 endif
260
261 if (con_spx(con_spx(je -1) +js).eq.bn) then
262 ib = i
263 jb = j
264 kb = k
265 endif
266 enddo
267 enddo
268 enddo
269
270 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
271 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
272
273 do l = 1,nn
274 is = nn*nn*(1 -1) +nn*(nn -1) +l
275 js = jsa +dljs*(l -1)
276 con_spx(con_spx(ie -1) +is) = &
277 con_spx(con_spx(je -1) +js)
278 enddo
279
280 else
281 do i = 2,(nn -1)
282 nnode = nnode +1
283 is = nn*nn*(1 -1) +nn*(nn -1) +i
284 con_spx(con_spx(ie -1) + is) = nnode
285 enddo
286 endif
287
288
289! Fourth edge down
290
291 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +1)
292 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +1)
293
294 je = ie
295 do i = node_pntr(an -1),node_pntr(an) -1
296 do j = node_pntr(bn -1),node_pntr(bn) -1
297 if (node_pntr(i).eq.node_pntr(j)) then
298 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
299 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
300 je = node_pntr(i)
301 endif
302 endif
303 enddo
304 enddo
305
306 if (je.ne.ie) then
307 ia = 0
308 ja = 0
309 ka = 0
310 ib = 0
311 jb = 0
312 kb = 0
313
314 do k = 1,nn
315 do j = 1,nn
316 do i = 1,nn
317 js = nn*nn*(k -1) +nn*(j -1) +i
318
319 if (con_spx(con_spx(je -1) +js).eq.an) then
320 ia = i
321 ja = j
322 ka = k
323 endif
324
325 if (con_spx(con_spx(je -1) +js).eq.bn) then
326 ib = i
327 jb = j
328 kb = k
329 endif
330 enddo
331 enddo
332 enddo
333
334 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
335 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
336
337 do l = 1,nn
338 is = nn*nn*(1 -1) +nn*(l -1) +1
339 js = jsa +dljs*(l -1)
340 con_spx(con_spx(ie -1) +is) = &
341 con_spx(con_spx(je -1) +js)
342 enddo
343
344 else
345 do j = 2,(nn -1)
346 nnode = nnode +1
347 is = nn*nn*(1 -1) +nn*(j -1) +1
348 con_spx(con_spx(ie -1) + is) = nnode
349 enddo
350 endif
351
352
353! First edge side
354
355 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +1)
356 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +1)
357
358 je = ie
359 do i = node_pntr(an -1),node_pntr(an) -1
360 do j = node_pntr(bn -1),node_pntr(bn) -1
361 if (node_pntr(i).eq.node_pntr(j)) then
362 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
363 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
364 je = node_pntr(i)
365 endif
366 endif
367 enddo
368 enddo
369
370 if (je.ne.ie) then
371 ia = 0
372 ja = 0
373 ka = 0
374 ib = 0
375 jb = 0
376 kb = 0
377
378 do k = 1,nn
379 do j = 1,nn
380 do i = 1,nn
381 js = nn*nn*(k -1) +nn*(j -1) +i
382
383 if (con_spx(con_spx(je -1) +js).eq.an) then
384 ia = i
385 ja = j
386 ka = k
387 endif
388
389 if (con_spx(con_spx(je -1) +js).eq.bn) then
390 ib = i
391 jb = j
392 kb = k
393 endif
394 enddo
395 enddo
396 enddo
397
398 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
399 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
400
401 do l = 1,nn
402 is = nn*nn*(l -1) +nn*(1 -1) +1
403 js = jsa +dljs*(l -1)
404 con_spx(con_spx(ie -1) +is) = &
405 con_spx(con_spx(je -1) +js)
406 enddo
407
408 else
409 do k = 2,(nn -1)
410 nnode = nnode +1
411 is = nn*nn*(k -1) +nn*(1 -1) +1
412 con_spx(con_spx(ie -1) + is) = nnode
413 enddo
414 endif
415
416
417! Second edge side
418
419 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +nn)
420 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +nn)
421
422 je = ie
423 do i = node_pntr(an -1),node_pntr(an) -1
424 do j = node_pntr(bn -1),node_pntr(bn) -1
425 if (node_pntr(i).eq.node_pntr(j)) then
426 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
427 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
428 je = node_pntr(i)
429 endif
430 endif
431 enddo
432 enddo
433
434 if (je.ne.ie) then
435 ia = 0
436 ja = 0
437 ka = 0
438 ib = 0
439 jb = 0
440 kb = 0
441
442 do k = 1,nn
443 do j = 1,nn
444 do i = 1,nn
445 js = nn*nn*(k -1) +nn*(j -1) +i
446
447 if (con_spx(con_spx(je -1) +js).eq.an) then
448 ia = i
449 ja = j
450 ka = k
451 endif
452
453 if (con_spx(con_spx(je -1) +js).eq.bn) then
454 ib = i
455 jb = j
456 kb = k
457 endif
458 enddo
459 enddo
460 enddo
461
462 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
463 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
464
465 do l = 1,nn
466 is = nn*nn*(l -1) +nn*(1 -1) +nn
467 js = jsa +dljs*(l -1)
468 con_spx(con_spx(ie -1) +is) = &
469 con_spx(con_spx(je -1) +js)
470 enddo
471
472 else
473 do k = 2,(nn -1)
474 nnode = nnode +1
475 is = nn*nn*(k -1) +nn*(1 -1) +nn
476 con_spx(con_spx(ie -1) + is) = nnode
477 enddo
478 endif
479
480
481! Third edge side
482
483 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +nn)
484 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +nn)
485
486 je = ie
487 do i = node_pntr(an -1),node_pntr(an) -1
488 do j = node_pntr(bn -1),node_pntr(bn) -1
489 if (node_pntr(i).eq.node_pntr(j)) then
490 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
491 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
492 je = node_pntr(i)
493 endif
494 endif
495 enddo
496 enddo
497
498 if (je.ne.ie) then
499 ia = 0
500 ja = 0
501 ka = 0
502 ib = 0
503 jb = 0
504 kb = 0
505
506 do k = 1,nn
507 do j = 1,nn
508 do i = 1,nn
509 js = nn*nn*(k -1) +nn*(j -1) +i
510
511 if (con_spx(con_spx(je -1) +js).eq.an) then
512 ia = i
513 ja = j
514 ka = k
515 endif
516
517 if (con_spx(con_spx(je -1) +js).eq.bn) then
518 ib = i
519 jb = j
520 kb = k
521 endif
522 enddo
523 enddo
524 enddo
525
526 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
527 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
528
529 do l = 1,nn
530 is = nn*nn*(l -1) +nn*(nn -1) +nn
531 js = jsa +dljs*(l -1)
532 con_spx(con_spx(ie -1) +is) = &
533 con_spx(con_spx(je -1) +js)
534 enddo
535
536 else
537 do k = 2,(nn -1)
538 nnode = nnode +1
539 is = nn*nn*(k -1) +nn*(nn -1) +nn
540 con_spx(con_spx(ie -1) + is) = nnode
541 enddo
542 endif
543
544
545! Fourth edge side
546
547 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +1)
548 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +1)
549
550 je = ie
551 do i = node_pntr(an -1),node_pntr(an) -1
552 do j = node_pntr(bn -1),node_pntr(bn) -1
553 if (node_pntr(i).eq.node_pntr(j)) then
554 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
555 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
556 je = node_pntr(i)
557 endif
558 endif
559 enddo
560 enddo
561
562 if (je.ne.ie) then
563 ia = 0
564 ja = 0
565 ka = 0
566
567 ib = 0
568 jb = 0
569 kb = 0
570
571 do k = 1,nn
572 do j = 1,nn
573 do i = 1,nn
574 js = nn*nn*(k -1) +nn*(j -1) +i
575
576 if (con_spx(con_spx(je -1) +js).eq.an) then
577 ia = i
578 ja = j
579 ka = k
580 endif
581
582 if (con_spx(con_spx(je -1) +js).eq.bn) then
583 ib = i
584 jb = j
585 kb = k
586 endif
587 enddo
588 enddo
589 enddo
590
591 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
592 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
593
594 do l = 1,nn
595 is = nn*nn*(l -1) +nn*(nn -1) +1
596 js = jsa +dljs*(l -1)
597 con_spx(con_spx(ie -1) +is) = &
598 con_spx(con_spx(je -1) +js)
599 enddo
600
601 else
602 do k = 2,(nn -1)
603 nnode = nnode +1
604 is = nn*nn*(k -1) +nn*(nn -1) +1
605 con_spx(con_spx(ie -1) + is) = nnode
606 enddo
607 endif
608
609
610! First edge up
611
612 an = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +1)
613 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +nn)
614
615 je = ie
616 do i = node_pntr(an -1),node_pntr(an) -1
617 do j = node_pntr(bn -1),node_pntr(bn) -1
618 if (node_pntr(i).eq.node_pntr(j)) then
619 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
620 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
621 je = node_pntr(i)
622 endif
623 endif
624 enddo
625 enddo
626
627 if (je.ne.ie) then
628 ia = 0
629 ja = 0
630 ka = 0
631
632 ib = 0
633 jb = 0
634 kb = 0
635
636 do k = 1,nn
637 do j = 1,nn
638 do i = 1,nn
639 js = nn*nn*(k -1) +nn*(j -1) +i
640 if (con_spx(con_spx(je -1) +js).eq.an) then
641 ia = i
642 ja = j
643 ka = k
644 endif
645
646 if (con_spx(con_spx(je -1) +js).eq.bn) then
647 ib = i
648 jb = j
649 kb = k
650 endif
651 enddo
652 enddo
653 enddo
654
655 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
656 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
657
658 do l = 1,nn
659 is = nn*nn*(nn -1) +nn*(1 -1) +l
660 js = jsa +dljs*(l -1)
661 con_spx(con_spx(ie -1) +is) = &
662 con_spx(con_spx(je -1) +js)
663 enddo
664
665 else
666 do i = 2,(nn -1)
667 nnode = nnode +1
668 is = nn*nn*(nn -1) +nn*(1 -1) +i
669 con_spx(con_spx(ie -1) + is) = nnode
670 enddo
671 endif
672
673
674! Second edge up
675
676 an = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +nn)
677 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +nn)
678
679 je = ie
680 do i = node_pntr(an -1),node_pntr(an) -1
681 do j = node_pntr(bn -1),node_pntr(bn) -1
682 if (node_pntr(i).eq.node_pntr(j)) then
683 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
684 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
685 je = node_pntr(i)
686 endif
687 endif
688 enddo
689 enddo
690
691 if (je.ne.ie) then
692 ia = 0
693 ja = 0
694 ka = 0
695 ib = 0
696 jb = 0
697 kb = 0
698
699 do k = 1,nn
700 do j = 1,nn
701 do i = 1,nn
702 js = nn*nn*(k -1) +nn*(j -1) +i
703
704 if (con_spx(con_spx(je -1) +js).eq.an) then
705 ia = i
706 ja = j
707 ka = k
708 endif
709
710 if (con_spx(con_spx(je -1) +js).eq.bn) then
711 ib = i
712 jb = j
713 kb = k
714 endif
715 enddo
716 enddo
717 enddo
718
719 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
720 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
721
722 do l = 1,nn
723 is = nn*nn*(nn -1) +nn*(l -1) +nn
724 js = jsa +dljs*(l -1)
725 con_spx(con_spx(ie -1) +is) = &
726 con_spx(con_spx(je -1) +js)
727 enddo
728
729 else
730 do j = 2,(nn -1)
731 nnode = nnode +1
732 is = nn*nn*(nn -1) +nn*(j -1) +nn
733 con_spx(con_spx(ie -1) + is) = nnode
734 enddo
735 endif
736
737
738! Third edge up
739
740 an = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +1)
741 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +nn)
742
743 je = ie
744 do i = node_pntr(an -1),node_pntr(an) -1
745 do j = node_pntr(bn -1),node_pntr(bn) -1
746 if (node_pntr(i).eq.node_pntr(j)) then
747 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
748 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
749 je = node_pntr(i)
750 endif
751 endif
752 enddo
753 enddo
754
755
756 if (je.ne.ie) then
757 ia = 0
758 ja = 0
759 ka = 0
760 ib = 0
761 jb = 0
762 kb = 0
763
764 do k = 1,nn
765 do j = 1,nn
766 do i = 1,nn
767 js = nn*nn*(k -1) +nn*(j -1) +i
768
769 if (con_spx(con_spx(je -1) +js).eq.an) then
770 ia = i
771 ja = j
772 ka = k
773 endif
774
775 if (con_spx(con_spx(je -1) +js).eq.bn) then
776 ib = i
777 jb = j
778 kb = k
779 endif
780 enddo
781 enddo
782 enddo
783
784 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
785 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
786
787 do l = 1,nn
788 is = nn*nn*(nn -1) +nn*(nn -1) +l
789 js = jsa +dljs*(l -1)
790 con_spx(con_spx(ie -1) +is) = &
791 con_spx(con_spx(je -1) +js)
792 enddo
793
794 else
795 do i = 2,(nn -1)
796 nnode = nnode +1
797 is = nn*nn*(nn -1) +nn*(nn -1) +i
798 con_spx(con_spx(ie -1) + is) = nnode
799 enddo
800 endif
801
802
803! Fourth edge up
804
805 an = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +1)
806 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +1)
807
808 je = ie
809 do i = node_pntr(an -1),node_pntr(an) -1
810 do j = node_pntr(bn -1),node_pntr(bn) -1
811 if (node_pntr(i).eq.node_pntr(j)) then
812 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
813 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
814 je = node_pntr(i)
815 endif
816 endif
817 enddo
818 enddo
819
820
821 if (je.ne.ie) then
822 ia = 0
823 ja = 0
824 ka = 0
825 ib = 0
826 jb = 0
827 kb = 0
828
829 do k = 1,nn
830 do j = 1,nn
831 do i = 1,nn
832 js = nn*nn*(k -1) +nn*(j -1) +i
833
834 if (con_spx(con_spx(je -1) +js).eq.an) then
835 ia = i
836 ja = j
837 ka = k
838 endif
839
840 if (con_spx(con_spx(je -1) +js).eq.bn) then
841 ib = i
842 jb = j
843 kb = k
844 endif
845 enddo
846 enddo
847 enddo
848
849 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
850 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
851
852 do l = 1,nn
853 is = nn*nn*(nn -1) +nn*(l -1) +1
854 js = jsa +dljs*(l -1)
855 con_spx(con_spx(ie -1) +is) = &
856 con_spx(con_spx(je -1) +js)
857 enddo
858
859 else
860 do j = 2,(nn -1)
861 nnode = nnode +1
862 is = nn*nn*(nn -1) +nn*(j -1) +1
863 con_spx(con_spx(ie -1) + is) = nnode
864 enddo
865 endif
866
867
868! Then construct face connectivity
869
870! Face down
871
872 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +1)
873 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +nn)
874 cn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +1)
875
876 je = ie
877 do i = node_pntr(an -1),node_pntr(an) -1
878 do j = node_pntr(bn -1),node_pntr(bn) -1
879 do k = node_pntr(cn -1),node_pntr(cn) -1
880 if ((node_pntr(i).eq.node_pntr(j))&
881 .and.(node_pntr(i).eq.node_pntr(k))) then
882 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
883 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
884 je = node_pntr(i)
885 endif
886 endif
887 enddo
888 enddo
889 enddo
890
891 if (je.ne.ie) then
892 ia = 0
893 ja = 0
894 ka = 0
895 ib = 0
896 jb = 0
897 kb = 0
898 ic = 0
899 jc = 0
900 kc = 0
901
902 do k = 1,nn
903 do j = 1,nn
904 do i = 1,nn
905 js = nn*nn*(k -1) +nn*(j -1) +i
906
907 if (con_spx(con_spx(je -1) +js).eq.an) then
908 ia = i
909 ja = j
910 ka = k
911 endif
912
913 if (con_spx(con_spx(je -1) +js).eq.bn) then
914 ib = i
915 jb = j
916 kb = k
917 endif
918
919 if (con_spx(con_spx(je -1) +js).eq.cn) then
920 ic = i
921 jc = j
922 kc = k
923 endif
924 enddo
925 enddo
926 enddo
927
928 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
929 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
930 dmjs = (nn*nn*(kc -ka) +nn*(jc -ja) +(ic -ia)) /(nn -1)
931
932 do m = 1,nn
933 do l = 1,nn
934 is = nn*nn*(1 -1) +nn*(m -1) +l
935 js = jsa +dmjs*(m -1) +dljs*(l -1)
936 con_spx(con_spx(ie -1) +is) = &
937 con_spx(con_spx(je -1) +js)
938 enddo
939 enddo
940
941 else
942 do m = 2,(nn -1)
943 do l = 2,(nn -1)
944 nnode = nnode +1
945 is = nn*nn*(1 -1) +nn*(m -1) +l
946 con_spx(con_spx(ie -1) + is) = nnode
947 enddo
948 enddo
949 endif
950
951
952! First face side
953
954 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +1)
955 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +nn)
956 cn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +1)
957
958 je = ie
959 do i = node_pntr(an -1),node_pntr(an) -1
960 do j = node_pntr(bn -1),node_pntr(bn) -1
961 do k = node_pntr(cn -1),node_pntr(cn) -1
962 if ((node_pntr(i).eq.node_pntr(j))&
963 .and.(node_pntr(i).eq.node_pntr(k))) then
964 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
965 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
966 je = node_pntr(i)
967 endif
968 endif
969 enddo
970 enddo
971 enddo
972
973 if (je.ne.ie) then
974 ia = 0
975 ja = 0
976 ka = 0
977 ib = 0
978 jb = 0
979 kb = 0
980 ic = 0
981 jc = 0
982 kc = 0
983
984 do k = 1,nn
985 do j = 1,nn
986 do i = 1,nn
987 js = nn*nn*(k -1) +nn*(j -1) +i
988
989 if (con_spx(con_spx(je -1) +js).eq.an) then
990 ia = i
991 ja = j
992 ka = k
993 endif
994
995 if (con_spx(con_spx(je -1) +js).eq.bn) then
996 ib = i
997 jb = j
998 kb = k
999 endif
1000
1001 if (con_spx(con_spx(je -1) +js).eq.cn) then
1002 ic = i
1003 jc = j
1004 kc = k
1005 endif
1006 enddo
1007 enddo
1008 enddo
1009
1010 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
1011 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
1012 dmjs = (nn*nn*(kc -ka) +nn*(jc -ja) +(ic -ia)) /(nn -1)
1013
1014 do m = 1,nn
1015 do l = 1,nn
1016 is = nn*nn*(m -1) +nn*(1 -1) +l
1017 js = jsa +dmjs*(m -1) +dljs*(l -1)
1018 con_spx(con_spx(ie -1) +is) = &
1019 con_spx(con_spx(je -1) +js)
1020 enddo
1021 enddo
1022
1023 else
1024 do m = 2,(nn -1)
1025 do l = 2,(nn -1)
1026 nnode = nnode +1
1027 is = nn*nn*(m -1) +nn*(1 -1) +l
1028 con_spx(con_spx(ie -1) + is) = nnode
1029 enddo
1030 enddo
1031 endif
1032
1033
1034! Second face side
1035
1036 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +nn)
1037 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +nn)
1038 cn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +nn)
1039
1040 je = ie
1041 do i = node_pntr(an -1),node_pntr(an) -1
1042 do j = node_pntr(bn -1),node_pntr(bn) -1
1043 do k = node_pntr(cn -1),node_pntr(cn) -1
1044 if ((node_pntr(i).eq.node_pntr(j))&
1045 .and.(node_pntr(i).eq.node_pntr(k))) then
1046 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
1047 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
1048 je = node_pntr(i)
1049 endif
1050 endif
1051 enddo
1052 enddo
1053 enddo
1054
1055 if (je.ne.ie) then
1056 ia = 0
1057 ja = 0
1058 ka = 0
1059 ib = 0
1060 jb = 0
1061 kb = 0
1062 ic = 0
1063 jc = 0
1064 kc = 0
1065
1066 do k = 1,nn
1067 do j = 1,nn
1068 do i = 1,nn
1069 js = nn*nn*(k -1) +nn*(j -1) +i
1070
1071 if (con_spx(con_spx(je -1) +js).eq.an) then
1072 ia = i
1073 ja = j
1074 ka = k
1075 endif
1076
1077 if (con_spx(con_spx(je -1) +js).eq.bn) then
1078 ib = i
1079 jb = j
1080 kb = k
1081 endif
1082
1083 if (con_spx(con_spx(je -1) +js).eq.cn) then
1084 ic = i
1085 jc = j
1086 kc = k
1087 endif
1088 enddo
1089 enddo
1090 enddo
1091
1092 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
1093 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
1094 dmjs = (nn*nn*(kc -ka) +nn*(jc -ja) +(ic -ia)) /(nn -1)
1095
1096 do m = 1,nn
1097 do l = 1,nn
1098 is = nn*nn*(m -1) +nn*(l -1) +nn
1099 js = jsa +dmjs*(m -1) +dljs*(l -1)
1100 con_spx(con_spx(ie -1) +is) = &
1101 con_spx(con_spx(je -1) +js)
1102 enddo
1103 enddo
1104
1105 else
1106 do m = 2,(nn -1)
1107 do l = 2,(nn -1)
1108 nnode = nnode +1
1109 is = nn*nn*(m -1) +nn*(l -1) +nn
1110 con_spx(con_spx(ie -1) + is) = nnode
1111 enddo
1112 enddo
1113 endif
1114
1115
1116! Third face side
1117
1118 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +1)
1119 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +nn)
1120 cn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +1)
1121
1122 je = ie
1123 do i = node_pntr(an -1),node_pntr(an) -1
1124 do j = node_pntr(bn -1),node_pntr(bn) -1
1125 do k = node_pntr(cn -1),node_pntr(cn) -1
1126 if ((node_pntr(i).eq.node_pntr(j))&
1127 .and.(node_pntr(i).eq.node_pntr(k))) then
1128 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
1129 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
1130 je = node_pntr(i)
1131 endif
1132 endif
1133 enddo
1134 enddo
1135 enddo
1136
1137 if (je.ne.ie) then
1138 ia = 0
1139 ja = 0
1140 ka = 0
1141 ib = 0
1142 jb = 0
1143 kb = 0
1144 ic = 0
1145 jc = 0
1146 kc = 0
1147
1148 do k = 1,nn
1149 do j = 1,nn
1150 do i = 1,nn
1151 js = nn*nn*(k -1) +nn*(j -1) +i
1152
1153 if (con_spx(con_spx(je -1) +js).eq.an) then
1154 ia = i
1155 ja = j
1156 ka = k
1157 endif
1158
1159 if (con_spx(con_spx(je -1) +js).eq.bn) then
1160 ib = i
1161 jb = j
1162 kb = k
1163 endif
1164
1165 if (con_spx(con_spx(je -1) +js).eq.cn) then
1166 ic = i
1167 jc = j
1168 kc = k
1169 endif
1170 enddo
1171 enddo
1172 enddo
1173
1174 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
1175 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
1176 dmjs = (nn*nn*(kc -ka) +nn*(jc -ja) +(ic -ia)) /(nn -1)
1177
1178 do m = 1,nn
1179 do l = 1,nn
1180 is = nn*nn*(m -1) +nn*(nn -1) +l
1181 js = jsa +dmjs*(m -1) +dljs*(l -1)
1182 con_spx(con_spx(ie -1) +is) = &
1183 con_spx(con_spx(je -1) +js)
1184 enddo
1185 enddo
1186
1187 else
1188 do m = 2,(nn -1)
1189 do l = 2,(nn -1)
1190 nnode = nnode +1
1191 is = nn*nn*(m -1) +nn*(nn -1) +l
1192 con_spx(con_spx(ie -1) + is) = nnode
1193 enddo
1194 enddo
1195 endif
1196
1197
1198! Fourth face side
1199
1200 an = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(1 -1) +1)
1201 bn = con_spx(con_spx(ie -1) +nn*nn*(1 -1) +nn*(nn -1) +1)
1202 cn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +1)
1203
1204 je = ie
1205 do i = node_pntr(an -1),node_pntr(an) -1
1206 do j = node_pntr(bn -1),node_pntr(bn) -1
1207 do k = node_pntr(cn -1),node_pntr(cn) -1
1208 if ((node_pntr(i).eq.node_pntr(j))&
1209 .and.(node_pntr(i).eq.node_pntr(k))) then
1210 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
1211 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
1212 je = node_pntr(i)
1213 endif
1214 endif
1215 enddo
1216 enddo
1217 enddo
1218
1219 if (je.ne.ie) then
1220 ia = 0
1221 ja = 0
1222 ka = 0
1223 ib = 0
1224 jb = 0
1225 kb = 0
1226 ic = 0
1227 jc = 0
1228 kc = 0
1229
1230 do k = 1,nn
1231 do j = 1,nn
1232 do i = 1,nn
1233 js = nn*nn*(k -1) +nn*(j -1) +i
1234
1235 if (con_spx(con_spx(je -1) +js).eq.an) then
1236 ia = i
1237 ja = j
1238 ka = k
1239 endif
1240
1241 if (con_spx(con_spx(je -1) +js).eq.bn) then
1242 ib = i
1243 jb = j
1244 kb = k
1245 endif
1246
1247 if (con_spx(con_spx(je -1) +js).eq.cn) then
1248 ic = i
1249 jc = j
1250 kc = k
1251 endif
1252 enddo
1253 enddo
1254 enddo
1255
1256 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
1257 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
1258 dmjs = (nn*nn*(kc -ka) +nn*(jc -ja) +(ic -ia)) /(nn -1)
1259
1260 do m = 1,nn
1261 do l = 1,nn
1262 is = nn*nn*(m -1) +nn*(l -1) +1
1263 js = jsa +dmjs*(m -1) +dljs*(l -1)
1264 con_spx(con_spx(ie -1) +is) = &
1265 con_spx(con_spx(je -1) +js)
1266 enddo
1267 enddo
1268
1269 else
1270 do m = 2,(nn -1)
1271 do l = 2,(nn -1)
1272 nnode = nnode +1
1273 is = nn*nn*(m -1) +nn*(l -1) +1
1274 con_spx(con_spx(ie -1) + is) = nnode
1275 enddo
1276 enddo
1277 endif
1278
1279
1280! Face up
1281
1282 an = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +1)
1283 bn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(1 -1) +nn)
1284 cn = con_spx(con_spx(ie -1) +nn*nn*(nn -1) +nn*(nn -1) +1)
1285
1286 je = ie
1287 do i = node_pntr(an -1),node_pntr(an) -1
1288 do j = node_pntr(bn -1),node_pntr(bn) -1
1289 do k = node_pntr(cn -1),node_pntr(cn) -1
1290 if ((node_pntr(i).eq.node_pntr(j))&
1291 .and.(node_pntr(i).eq.node_pntr(k))) then
1292 if ((node_pntr(i).lt.ie).and.((nn*nn*nn +1)&
1293 .eq.(con_spx(node_pntr(i)) - con_spx(node_pntr(i) -1)))) then
1294 je = node_pntr(i)
1295 endif
1296 endif
1297 enddo
1298 enddo
1299 enddo
1300
1301 if (je.ne.ie) then
1302 ia = 0
1303 ja = 0
1304 ka = 0
1305 ib = 0
1306 jb = 0
1307 kb = 0
1308 ic = 0
1309 jc = 0
1310 kc = 0
1311
1312 do k = 1,nn
1313 do j = 1,nn
1314 do i = 1,nn
1315 js = nn*nn*(k -1) +nn*(j -1) +i
1316
1317 if (con_spx(con_spx(je -1) +js).eq.an) then
1318 ia = i
1319 ja = j
1320 ka = k
1321 endif
1322
1323 if (con_spx(con_spx(je -1) +js).eq.bn) then
1324 ib = i
1325 jb = j
1326 kb = k
1327 endif
1328
1329 if (con_spx(con_spx(je -1) +js).eq.cn) then
1330 ic = i
1331 jc = j
1332 kc = k
1333 endif
1334 enddo
1335 enddo
1336 enddo
1337
1338 jsa = nn*nn*(ka -1) +nn*(ja -1) +(ia -1) +1
1339 dljs = (nn*nn*(kb -ka) +nn*(jb -ja) +(ib -ia)) /(nn -1)
1340 dmjs = (nn*nn*(kc -ka) +nn*(jc -ja) +(ic -ia)) /(nn -1)
1341
1342 do m = 1,nn
1343 do l = 1,nn
1344 is = nn*nn*(nn -1) +nn*(m -1) +l
1345 js = jsa +dmjs*(m -1) +dljs*(l -1)
1346 con_spx(con_spx(ie -1) +is) = &
1347 con_spx(con_spx(je -1) +js)
1348 enddo
1349 enddo
1350
1351 else
1352 do m = 2,(nn -1)
1353 do l = 2,(nn -1)
1354 nnode = nnode +1
1355 is = nn*nn*(nn -1) +nn*(m -1) +l
1356 con_spx(con_spx(ie -1) + is) = nnode
1357 enddo
1358 enddo
1359 endif
1360
1361
1362! End of edge and face connectivity
1363
1364! Finally fill inside the elements
1365
1366 do k = 2,(nn -1)
1367 do j = 2,(nn -1)
1368 do i = 2,(nn -1)
1369 nnode = nnode +1
1370 is = nn*nn*(k -1) +nn*(j -1) +i
1371 con_spx(con_spx(ie -1) + is) = nnode
1372 enddo
1373 enddo
1374 enddo
1375
1376 endif
1377 enddo
1378 enddo
1379
1380 return
1381
1382 end subroutine make_spx_con
1383
make_spx_con
subroutine make_spx_con(nelem, con_mac, nmat, tag_mat, sdeg, nnz_pntr, node_pntr, con_nnz, c
Makes spectral connectivity vector.
Definition MAKE_SPX_CON.f90:34