#include "InfoStream.h"#include "ComputeCrossterms.h"#include "Molecule.h"#include "Parameters.h"#include "Node.h"#include "ReductionMgr.h"#include "Lattice.h"#include "PressureProfile.h"#include "Debug.h"Go to the source code of this file.
Defines | |
| #define | FASTER |
| #define | INDEX(ncols, i, j) ((i)*ncols + (j)) |
Enumerations | |
| enum | { CMAP_TABLE_DIM = 25, CMAP_SETUP_DIM = 24, CMAP_SPACING = 15, CMAP_PHI_0 = -180, CMAP_PSI_0 = -180 } |
Functions | |
| void | lu_decomp_nopivot (double *m, int n) |
| void | forward_back_sub (double *b, double *m, int n) |
| void | crossterm_setup (CrosstermData *table) |
|
|
Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved. Definition at line 17 of file ComputeCrossterms.C. |
|
|
Definition at line 27 of file ComputeCrossterms.C. Referenced by CrosstermElem::computeForce(), crossterm_setup(), forward_back_sub(), and lu_decomp_nopivot(). |
|
|
Definition at line 19 of file ComputeCrossterms.C. 00019 {
00020 CMAP_TABLE_DIM = 25,
00021 CMAP_SETUP_DIM = 24,
00022 CMAP_SPACING = 15,
00023 CMAP_PHI_0 = -180,
00024 CMAP_PSI_0 = -180
00025 };
|
|
|
Definition at line 539 of file ComputeCrossterms.C. References CMAP_SETUP_DIM, CMAP_SPACING, CMAP_TABLE_DIM, CrosstermData::d00, CrosstermData::d01, CrosstermData::d10, CrosstermData::d11, forward_back_sub(), INDEX, lu_decomp_nopivot(), and PI. 00540 {
00541 enum { D = CMAP_TABLE_DIM };
00542 enum { N = CMAP_SETUP_DIM };
00543 const double h_1 = 1.0 / ( CMAP_SPACING * PI / 180.0) ;
00544 const double h_2 = h_1 * h_1;
00545 const double tr_h = 3.0 * h_1;
00546 int i, j;
00547 int ij;
00548 int ijp1p1, ijm1m1, ijp1m1, ijm1p1;
00549 int ijp2p2, ijm2m2, ijp2m2, ijm2p2;
00550
00551 /* allocate spline coefficient matrix */
00552 double* const m = new double[N*N];
00553 memset(m,0,N*N*sizeof(double));
00554
00555 /* initialize spline coefficient matrix */
00556 m[0] = 4;
00557 for (i = 1; i < N; i++) {
00558 m[INDEX(N,i-1,i)] = 1;
00559 m[INDEX(N,i,i-1)] = 1;
00560 m[INDEX(N,i,i)] = 4;
00561 }
00562 m[INDEX(N,0,N-1)] = 1;
00563 m[INDEX(N,N-1,0)] = 1;
00564
00565 /* compute LU-decomposition for this matrix */
00566 lu_decomp_nopivot(m, N);
00567
00568 /* allocate vector for solving spline derivatives */
00569 double* const v = new double[N];
00570 memset(v,0,N*sizeof(double));
00571
00572 /* march through rows of table */
00573 for (i = 0; i < N; i++) {
00574
00575 /* setup RHS vector for solving spline derivatives */
00576 v[0] = tr_h * (table[INDEX(D,i,1)].d00 - table[INDEX(D,i,N-1)].d00);
00577 for (j = 1; j < N; j++) {
00578 v[j] = tr_h * (table[INDEX(D,i,j+1)].d00 - table[INDEX(D,i,j-1)].d00);
00579 }
00580
00581 /* solve system, returned into vector */
00582 forward_back_sub(v, m, N);
00583
00584 /* store values as derivatives wrt differenced table values */
00585 for (j = 0; j < N; j++) {
00586 table[INDEX(D,i,j)].d01 = v[j];
00587 }
00588 table[INDEX(D,i,N)].d01 = v[0];
00589 }
00590 for (j = 0; j <= N; j++) {
00591 table[INDEX(D,N,j)].d01 = table[INDEX(D,0,j)].d01;
00592 }
00593
00594 /* march through columns of table */
00595 for (j = 0; j < N; j++) {
00596
00597 /* setup RHS vector for solving spline derivatives */
00598 v[0] = tr_h * (table[INDEX(D,1,j)].d00 - table[INDEX(D,N-1,j)].d00);
00599 for (i = 1; i < N; i++) {
00600 v[i] = tr_h * (table[INDEX(D,i+1,j)].d00 - table[INDEX(D,i-1,j)].d00);
00601 }
00602
00603 /* solve system, returned into vector */
00604 forward_back_sub(v, m, N);
00605
00606 /* store values as derivatives wrt differenced table values */
00607 for (i = 0; i < N; i++) {
00608 table[INDEX(D,i,j)].d10 = v[i];
00609 }
00610 table[INDEX(D,N,j)].d10 = v[0];
00611 }
00612 for (i = 0; i <= N; i++) {
00613 table[INDEX(D,i,N)].d10 = table[INDEX(D,i,0)].d10;
00614 }
00615
00616 /* use 4th order difference approx for cross-derivative */
00617 for (i = 0; i < N; i++) {
00618 for (j = 0; j < N; j++) {
00619 ij = INDEX(D,i,j);
00620 ijp1p1 = INDEX(D,(i+1)%N,(j+1)%N);
00621 ijm1m1 = INDEX(D,(i+N-1)%N,(j+N-1)%N);
00622 ijp1m1 = INDEX(D,(i+1)%N,(j+N-1)%N);
00623 ijm1p1 = INDEX(D,(i+N-1)%N,(j+1)%N);
00624 ijp2p2 = INDEX(D,(i+2)%N,(j+2)%N);
00625 ijm2m2 = INDEX(D,(i+N-2)%N,(j+N-2)%N);
00626 ijp2m2 = INDEX(D,(i+2)%N,(j+N-2)%N);
00627 ijm2p2 = INDEX(D,(i+N-2)%N,(j+2)%N);
00628 table[ij].d11 = h_2 * (1./3) * (table[ijp1p1].d00 + table[ijm1m1].d00
00629 - table[ijp1m1].d00 - table[ijm1p1].d00
00630 - (1./16) * (table[ijp2p2].d00 + table[ijm2m2].d00
00631 - table[ijp2m2].d00 - table[ijm2p2].d00));
00632 }
00633 }
00634
00635 #if 0
00636 /* use 2nd order difference approx for cross-derivative */
00637 for (i = 0; i < N; i++) {
00638 for (j = 0; j < N; j++) {
00639 ij = INDEX(D,i,j);
00640 ijp1p1 = INDEX(D,(i+1)%N,(j+1)%N);
00641 ijm1m1 = INDEX(D,(i+N-1)%N,(j+N-1)%N);
00642 ijp1m1 = INDEX(D,(i+1)%N,(j+N-1)%N);
00643 ijm1p1 = INDEX(D,(i+N-1)%N,(j+1)%N);
00644 table[ij].d11 = h_2 * (1./4) * (table[ijp1p1].d00 + table[ijm1m1].d00
00645 - table[ijp1m1].d00 - table[ijm1p1].d00);
00646 }
00647 }
00648 #endif
00649
00650 /* fill in redundant edge values */
00651 for (i = 0; i < N; i++) {
00652 table[INDEX(D,i,N)].d11 = table[INDEX(D,i,0)].d11;
00653 table[INDEX(D,N,i)].d11 = table[INDEX(D,0,i)].d11;
00654 }
00655 table[INDEX(D,N,N)].d11 = table[INDEX(D,0,0)].d11;
00656
00657 /* done with temp storage */
00658 delete [] m;
00659 delete [] v;
00660
00661 }
|
|
||||||||||||||||
|
Definition at line 681 of file ComputeCrossterms.C. References INDEX. Referenced by crossterm_setup(). 00682 {
00683 int i, j;
00684
00685 /* in place forward elimination (solves Ly=b) using lower triangle of m */
00686 for (j = 0; j < n-1; j++) {
00687 for (i = j+1; i < n; i++) {
00688 b[i] -= m[INDEX(n,i,j)] * b[j];
00689 }
00690 }
00691 /* in place back substitution (solves Ux=y) using upper triangle of m */
00692 for (j = n-1; j >= 0; j--) {
00693 b[j] /= m[INDEX(n,j,j)];
00694 for (i = j-1; i >= 0; i--) {
00695 b[i] -= m[INDEX(n,i,j)] * b[j];
00696 }
00697 }
00698 }
|
|
||||||||||||
|
Definition at line 664 of file ComputeCrossterms.C. References INDEX. Referenced by crossterm_setup(). 00665 {
00666 double l_ik;
00667 int i, j, k;
00668
00669 for (k = 0; k < n-1; k++) {
00670 for (i = k+1; i < n; i++) {
00671 l_ik = m[INDEX(n,i,k)] / m[INDEX(n,k,k)];
00672 for (j = k; j < n; j++) {
00673 m[INDEX(n,i,j)] -= l_ik * m[INDEX(n,k,j)];
00674 }
00675 m[INDEX(n,i,k)] = l_ik;
00676 }
00677 }
00678 }
|
1.3.9.1