#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(), lu_decomp_nopivot(), Parameters::read_energy_type_bothcubspline(), and Parameters::read_energy_type_cubspline(). |
|
|
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 526 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, j, lu_decomp_nopivot(), and PI. 00527 {
00528 enum { D = CMAP_TABLE_DIM };
00529 enum { N = CMAP_SETUP_DIM };
00530 const double h_1 = 1.0 / ( CMAP_SPACING * PI / 180.0) ;
00531 const double h_2 = h_1 * h_1;
00532 const double tr_h = 3.0 * h_1;
00533 int i, j;
00534 int ij;
00535 int ijp1p1, ijm1m1, ijp1m1, ijm1p1;
00536 int ijp2p2, ijm2m2, ijp2m2, ijm2p2;
00537
00538 /* allocate spline coefficient matrix */
00539 double* const m = new double[N*N];
00540 memset(m,0,N*N*sizeof(double));
00541
00542 /* initialize spline coefficient matrix */
00543 m[0] = 4;
00544 for (i = 1; i < N; i++) {
00545 m[INDEX(N,i-1,i)] = 1;
00546 m[INDEX(N,i,i-1)] = 1;
00547 m[INDEX(N,i,i)] = 4;
00548 }
00549 m[INDEX(N,0,N-1)] = 1;
00550 m[INDEX(N,N-1,0)] = 1;
00551
00552 /* compute LU-decomposition for this matrix */
00553 lu_decomp_nopivot(m, N);
00554
00555 /* allocate vector for solving spline derivatives */
00556 double* const v = new double[N];
00557 memset(v,0,N*sizeof(double));
00558
00559 /* march through rows of table */
00560 for (i = 0; i < N; i++) {
00561
00562 /* setup RHS vector for solving spline derivatives */
00563 v[0] = tr_h * (table[INDEX(D,i,1)].d00 - table[INDEX(D,i,N-1)].d00);
00564 for (j = 1; j < N; j++) {
00565 v[j] = tr_h * (table[INDEX(D,i,j+1)].d00 - table[INDEX(D,i,j-1)].d00);
00566 }
00567
00568 /* solve system, returned into vector */
00569 forward_back_sub(v, m, N);
00570
00571 /* store values as derivatives wrt differenced table values */
00572 for (j = 0; j < N; j++) {
00573 table[INDEX(D,i,j)].d01 = v[j];
00574 }
00575 table[INDEX(D,i,N)].d01 = v[0];
00576 }
00577 for (j = 0; j <= N; j++) {
00578 table[INDEX(D,N,j)].d01 = table[INDEX(D,0,j)].d01;
00579 }
00580
00581 /* march through columns of table */
00582 for (j = 0; j < N; j++) {
00583
00584 /* setup RHS vector for solving spline derivatives */
00585 v[0] = tr_h * (table[INDEX(D,1,j)].d00 - table[INDEX(D,N-1,j)].d00);
00586 for (i = 1; i < N; i++) {
00587 v[i] = tr_h * (table[INDEX(D,i+1,j)].d00 - table[INDEX(D,i-1,j)].d00);
00588 }
00589
00590 /* solve system, returned into vector */
00591 forward_back_sub(v, m, N);
00592
00593 /* store values as derivatives wrt differenced table values */
00594 for (i = 0; i < N; i++) {
00595 table[INDEX(D,i,j)].d10 = v[i];
00596 }
00597 table[INDEX(D,N,j)].d10 = v[0];
00598 }
00599 for (i = 0; i <= N; i++) {
00600 table[INDEX(D,i,N)].d10 = table[INDEX(D,i,0)].d10;
00601 }
00602
00603 /* use 4th order difference approx for cross-derivative */
00604 for (i = 0; i < N; i++) {
00605 for (j = 0; j < N; j++) {
00606 ij = INDEX(D,i,j);
00607 ijp1p1 = INDEX(D,(i+1)%N,(j+1)%N);
00608 ijm1m1 = INDEX(D,(i+N-1)%N,(j+N-1)%N);
00609 ijp1m1 = INDEX(D,(i+1)%N,(j+N-1)%N);
00610 ijm1p1 = INDEX(D,(i+N-1)%N,(j+1)%N);
00611 ijp2p2 = INDEX(D,(i+2)%N,(j+2)%N);
00612 ijm2m2 = INDEX(D,(i+N-2)%N,(j+N-2)%N);
00613 ijp2m2 = INDEX(D,(i+2)%N,(j+N-2)%N);
00614 ijm2p2 = INDEX(D,(i+N-2)%N,(j+2)%N);
00615 table[ij].d11 = h_2 * (1./3) * (table[ijp1p1].d00 + table[ijm1m1].d00
00616 - table[ijp1m1].d00 - table[ijm1p1].d00
00617 - (1./16) * (table[ijp2p2].d00 + table[ijm2m2].d00
00618 - table[ijp2m2].d00 - table[ijm2p2].d00));
00619 }
00620 }
00621
00622 #if 0
00623 /* use 2nd order difference approx for cross-derivative */
00624 for (i = 0; i < N; i++) {
00625 for (j = 0; j < N; j++) {
00626 ij = INDEX(D,i,j);
00627 ijp1p1 = INDEX(D,(i+1)%N,(j+1)%N);
00628 ijm1m1 = INDEX(D,(i+N-1)%N,(j+N-1)%N);
00629 ijp1m1 = INDEX(D,(i+1)%N,(j+N-1)%N);
00630 ijm1p1 = INDEX(D,(i+N-1)%N,(j+1)%N);
00631 table[ij].d11 = h_2 * (1./4) * (table[ijp1p1].d00 + table[ijm1m1].d00
00632 - table[ijp1m1].d00 - table[ijm1p1].d00);
00633 }
00634 }
00635 #endif
00636
00637 /* fill in redundant edge values */
00638 for (i = 0; i < N; i++) {
00639 table[INDEX(D,i,N)].d11 = table[INDEX(D,i,0)].d11;
00640 table[INDEX(D,N,i)].d11 = table[INDEX(D,0,i)].d11;
00641 }
00642 table[INDEX(D,N,N)].d11 = table[INDEX(D,0,0)].d11;
00643
00644 /* done with temp storage */
00645 delete [] m;
00646 delete [] v;
00647
00648 }
|
|
||||||||||||||||
|
Definition at line 668 of file ComputeCrossterms.C. Referenced by crossterm_setup(). 00669 {
00670 int i, j;
00671
00672 /* in place forward elimination (solves Ly=b) using lower triangle of m */
00673 for (j = 0; j < n-1; j++) {
00674 for (i = j+1; i < n; i++) {
00675 b[i] -= m[INDEX(n,i,j)] * b[j];
00676 }
00677 }
00678 /* in place back substitution (solves Ux=y) using upper triangle of m */
00679 for (j = n-1; j >= 0; j--) {
00680 b[j] /= m[INDEX(n,j,j)];
00681 for (i = j-1; i >= 0; i--) {
00682 b[i] -= m[INDEX(n,i,j)] * b[j];
00683 }
00684 }
00685 }
|
|
||||||||||||
|
Definition at line 651 of file ComputeCrossterms.C. Referenced by crossterm_setup(). 00652 {
00653 double l_ik;
00654 int i, j, k;
00655
00656 for (k = 0; k < n-1; k++) {
00657 for (i = k+1; i < n; i++) {
00658 l_ik = m[INDEX(n,i,k)] / m[INDEX(n,k,k)];
00659 for (j = k; j < n; j++) {
00660 m[INDEX(n,i,j)] -= l_ik * m[INDEX(n,k,j)];
00661 }
00662 m[INDEX(n,i,k)] = l_ik;
00663 }
00664 }
00665 }
|
1.3.9.1