|
00001 // Copyright (C) 2004 Johan Hoffman and Anders Logg. 00002 // Licensed under the GNU GPL Version 2. 00003 00004 #ifndef __POISSONSYSTEM_H 00005 #define __POISSONSYSTEM_H 00006 00007 #include <dolfin/NewPDE.h> 00008 00009 namespace dolfin 00010 { 00011 00013 00014 class PoissonSystem : public NewPDE 00015 { 00016 public: 00017 00018 PoissonSystem(Function& f) : NewPDE(8, true, false), w0(8) 00019 { 00020 // Add functions 00021 add(w0, f); 00022 00023 // Set nonzero pattern (add the rest of the entries) 00024 nonzero.push_back(IndexPair(0,0)); 00025 } 00026 00027 unsigned int dim() const 00028 { 00029 return 3; 00030 } 00031 00032 unsigned int dof(unsigned int i, const Cell& cell) const 00033 { 00034 // Something needs to be done here for systems 00035 00036 return cell.nodeID(i); 00037 } 00038 00039 void interiorElementMatrix(NewArray<NewArray<real> >& A) const 00040 { 00041 real tmp0 = det / 6.0; 00042 00043 real G00 = tmp0*(g00*g00 + g01*g01 + g02*g02); 00044 real G01 = tmp0*(g00*g10 + g01*g11 + g02*g12); 00045 real G02 = tmp0*(g00*g20 + g01*g21 + g02*g22); 00046 real G11 = tmp0*(g10*g10 + g11*g11 + g12*g12); 00047 real G12 = tmp0*(g10*g20 + g11*g21 + g12*g22); 00048 real G22 = tmp0*(g20*g20 + g21*g21 + g22*g22); 00049 00050 A[1][1] = G00; 00051 A[1][2] = G01; 00052 A[1][3] = G02; 00053 A[2][2] = G11; 00054 A[2][3] = G12; 00055 A[3][3] = G22; 00056 A[0][1] = - A[1][1] - A[1][2] - A[1][3]; 00057 A[0][2] = - A[1][2] - A[2][2] - A[2][3]; 00058 A[0][3] = - A[1][3] - A[2][3] - A[3][3]; 00059 A[0][0] = - A[0][1] - A[0][2] - A[0][3]; 00060 A[1][0] = A[0][1]; 00061 A[2][0] = A[0][2]; 00062 A[2][1] = A[1][2]; 00063 A[3][0] = A[0][3]; 00064 A[3][1] = A[1][3]; 00065 A[3][2] = A[2][3]; 00066 00067 A[4][4] = A[0][0]; 00068 A[4][5] = A[0][1]; 00069 A[4][6] = A[0][2]; 00070 A[4][7] = A[0][3]; 00071 A[5][4] = A[1][0]; 00072 A[5][5] = A[1][1]; 00073 A[5][6] = A[1][2]; 00074 A[5][7] = A[1][3]; 00075 A[6][4] = A[2][0]; 00076 A[6][5] = A[2][1]; 00077 A[6][6] = A[2][2]; 00078 A[6][7] = A[2][3]; 00079 A[7][4] = A[3][0]; 00080 A[7][5] = A[3][1]; 00081 A[7][6] = A[3][2]; 00082 A[7][7] = A[3][3]; 00083 } 00084 00085 void interiorElementVector(NewArray<real>& b) const 00086 { 00087 real tmp0 = det / 120.0; 00088 00089 real G0 = tmp0*w0[0]; 00090 real G1 = tmp0*w0[1]; 00091 real G2 = tmp0*w0[2]; 00092 real G3 = tmp0*w0[3]; 00093 real G4 = tmp0*w0[0]; 00094 real G5 = tmp0*w0[1]; 00095 real G6 = tmp0*w0[2]; 00096 real G7 = tmp0*w0[3]; 00097 00098 real tmp1 = G0 + G1 + G2 + G3; 00099 00100 b[0] = tmp1 + G0; 00101 b[1] = tmp1 + G1; 00102 b[2] = tmp1 + G2; 00103 b[3] = tmp1 + G3; 00104 00105 real tmp2 = G4 + G5 + G6 + G7; 00106 00107 b[4] = tmp2 + G4; 00108 b[5] = tmp2 + G5; 00109 b[6] = tmp2 + G6; 00110 b[7] = tmp2 + G7; 00111 } 00112 00113 private: 00114 00115 NewArray<real> w0; 00116 00117 }; 00118 00119 } 00120 00121 #endif
Documentation automatically generated with Doxygen on 10 Sep 2004.