00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "SundanceVectorCalculus.hpp"
00032 #include "SundanceDerivative.hpp"
00033 #include "SundanceListExpr.hpp"
00034
00035
00036 namespace Sundance
00037 {
00038
00039 Expr gradient(int dim)
00040 {
00041 Array<Expr> rtn(dim);
00042 for (int i=0; i<dim; i++)
00043 {
00044 rtn[i] = new Derivative(i);
00045 }
00046 return new ListExpr(rtn);
00047 }
00048
00049 Expr div(const Expr& f)
00050 {
00051 Expr rtn = 0.0;
00052 for (int i=0; i<f.size(); i++)
00053 {
00054 Expr d = new Derivative(i);
00055 rtn = rtn + d*f[i];
00056 }
00057 return rtn;
00058 }
00059
00060
00061 Expr cross(const Expr& a, const Expr& b)
00062 {
00063 TEUCHOS_TEST_FOR_EXCEPTION(a.size() != b.size(), std::runtime_error,
00064 "mismatched vector sizes in cross(a,b): a.size()=" << a.size()
00065 << ", b.size()=" << b.size());
00066
00067 TEUCHOS_TEST_FOR_EXCEPTION(a.size() < 2 || a.size() > 3, std::runtime_error,
00068 "cross(a,b) undefined for dim=" << a.size());
00069
00070 if (a.size()==2)
00071 {
00072 return a[0]*b[1] - a[1]*b[0];
00073 }
00074 else
00075 {
00076 return List(
00077 cross(List(a[1],a[2]), List(b[1],b[2])),
00078 -cross(List(a[0],a[2]), List(b[0],b[2])),
00079 cross(List(a[0],a[1]), List(b[0],b[1]))
00080 );
00081 }
00082 }
00083
00084 Expr curl(const Expr& f)
00085 {
00086 Expr del = gradient(f.size());
00087
00088 return cross(del, f);
00089 }
00090
00091 Expr colonProduct(const Expr& A, const Expr& B)
00092 {
00093 int nA = 0;
00094 int nB = 0;
00095 TEUCHOS_TEST_FOR_EXCEPTION(!isSquareMatrix(A, nA), std::runtime_error,
00096 "Colon product expected argument A=" << A << " to be a square matrix");
00097 TEUCHOS_TEST_FOR_EXCEPTION(!isSquareMatrix(B, nB), std::runtime_error,
00098 "Colon product expected argument B=" << B << " to be a square matrix");
00099
00100 TEUCHOS_TEST_FOR_EXCEPTION(nA!=nB, std::runtime_error,
00101 "Colon product expected operands A=" << A << " and B=" << B
00102 << " to have identical sizes");
00103
00104 Expr rtn = 0.0;
00105 for (int i=0; i<nA; i++)
00106 {
00107 for (int j=0; j<nA; j++)
00108 {
00109 rtn = rtn + A[i][j] * B[i][j];
00110 }
00111 }
00112
00113 return rtn;
00114 }
00115
00116 Expr outerProduct(const Expr& A, const Expr& B)
00117 {
00118 int nA = 0;
00119 int nB = 0;
00120 TEUCHOS_TEST_FOR_EXCEPTION(!isVector(A, nA), std::runtime_error,
00121 "Outer product expected argument A=" << A << " to be a vector");
00122 TEUCHOS_TEST_FOR_EXCEPTION(!isVector(B, nB), std::runtime_error,
00123 "Outer product expected argument B=" << B << " to be a vector");
00124
00125 TEUCHOS_TEST_FOR_EXCEPTION(nA!=nB, std::runtime_error,
00126 "Colon product expected operands A=" << A << " and B=" << B
00127 << " to have identical sizes");
00128
00129 Array<Expr> rtn(nA);
00130 for (int i=0; i<nA; i++)
00131 {
00132 rtn[i] = A[i] * B;
00133 }
00134
00135 return new ListExpr(rtn);
00136 }
00137
00138 bool isVector(const Expr& x, int& N)
00139 {
00140 N = 0;
00141 if (x.size() == x.totalSize())
00142 {
00143 N = x.size();
00144 return true;
00145 }
00146 else
00147 {
00148 return false;
00149 }
00150 }
00151
00152
00153 bool isSquareMatrix(const Expr& x, int& N)
00154 {
00155 N = 0;
00156
00157 if (x.size()*x.size() != x.totalSize()) return false;
00158 N = x.size();
00159 for (int i=0; i<N; i++)
00160 {
00161 int M = 0;
00162 if (!isVector(x[i],M)) return false;
00163 if (M != N) return false;
00164 }
00165 return true;
00166 }
00167
00168
00169
00170
00171 }
00172