AFEPack
|
00001 00011 #ifndef __BoundaryCondition_h__ 00012 #define __BoundaryCondition_h__ 00013 00014 #include <cstdarg> 00015 #include <vector> 00016 #include <set> 00017 #include <map> 00018 00019 #include <lac/sparse_matrix.h> 00020 #include <lac/vector.h> 00021 00022 #include "Miscellaneous.h" 00023 00037 class BCondition { 00038 public: 00040 static const int DIRICHLET; 00041 static const int NEUMANN; 00042 static const int ROBIN; 00043 00044 public: 00045 /*@\{*/ 00047 BCondition() {} 00048 BCondition(int type) : _type(type) {} 00049 BCondition(const BCondition& b) : 00050 _type(b._type), _bmark(b._bmark) {} 00051 virtual ~BCondition() {} 00052 /*@\}*/ 00053 00054 public: 00056 BCondition& operator=(const BCondition& b) { 00057 _type = b._type; 00058 _bmark = b._bmark; 00059 00060 return *this; 00061 } 00062 00063 /*@\{*/ 00065 int type() const {return _type;} 00066 int& type() {return _type;} 00067 /*@\}*/ 00068 00069 /*@\{*/ 00071 const std::set<int, std::less<int> >& bound_mark() const {return _bmark;} 00072 std::set<int, std::less<int> >& bound_mark() {return _bmark;} 00073 /*@\}*/ 00074 00075 /*@\{*/ 00077 void add_mark(int bm) {add_one_mark(bm);} 00078 void add_mark(const std::vector<int>& bm) { 00079 u_int n = bm.size(); 00080 for (u_int i = 0;i < n;++ i) { 00081 add_one_mark(bm[i]); 00082 } 00083 } 00084 template <class V> 00085 void add_mark(u_int n, const V& bm) { 00086 for (u_int i = 0;i < n;i ++) { 00087 add_one_mark(bm[i]); 00088 } 00089 } 00090 void add_mark(u_int n, int bm0, int bm1, ...) { 00091 add_one_mark(bm0); 00092 add_one_mark(bm1); 00093 00094 va_list ap; 00095 va_start(ap, bm1); 00096 for (u_int i = 2;i < n;i ++) { 00097 add_one_mark(va_arg(ap, int)); 00098 } 00099 va_end(ap); 00100 } 00101 /*@\}*/ 00102 00103 private: 00104 void add_one_mark(int bm) { 00105 if (bm == 0) { 00106 std::cerr << "棺ϱʶ 0 ʾڲܼӸ߽" 00107 << std::endl; 00108 return; 00109 } 00110 _bmark.insert(bm); 00111 } 00112 00113 public: 00115 virtual void value(const void * p, void * v) const {} 00117 virtual void gradient(const void * p, void * g) const {} 00118 00119 private: 00120 int _type; 00121 std::set<int, std::less<int> > _bmark; 00122 }; 00123 00124 00125 00130 template <class P, class V, class G = std::vector<V> > 00131 class BCFunction : public BCondition { 00132 public: 00133 typedef void (*val_fun_t)(const P&, V&); 00134 typedef void (*grad_fun_t)(const P&, G&); 00135 00136 private: 00137 static void _default_val(const P& p, V& v) {v = V();} 00138 static void _default_grad(const P& p, G& g) {g = G();} 00139 00140 val_fun_t _val_fun; 00141 grad_fun_t _grad_fun; 00142 00143 public: 00144 BCFunction(val_fun_t vf = &_default_val, 00145 grad_fun_t gf = &_default_grad) : 00146 _val_fun(vf), 00147 _grad_fun(gf) {} 00148 00149 BCFunction(int type, 00150 val_fun_t vf = &_default_val, 00151 grad_fun_t gf = &_default_grad) : 00152 BCondition(type), 00153 _val_fun(vf), 00154 _grad_fun(gf) {} 00155 00156 virtual ~BCFunction() {} 00157 00158 public: 00159 /*@\{*/ 00161 const val_fun_t& value_fun_ptr() const { 00162 return _val_fun; 00163 } 00164 val_fun_t& value_fun_ptr() { 00165 return _val_fun; 00166 } 00167 00168 const grad_fun_t& gradient_fun_ptr() const { 00169 return _grad_fun; 00170 } 00171 grad_fun_t& gradient_fun_ptr() { 00172 return _grad_fun; 00173 } 00174 /*@\}*/ 00175 00176 public: 00177 /*@\{*/ 00179 virtual void value(const void * p, void * v) const { 00180 (*_val_fun)(*(const P *)p, *(V *)v); 00181 } 00182 virtual void gradient(const void * p, void * g) const { 00183 (*_grad_fun)(*(const P *)p, *(G *)g); 00184 } 00185 /*@\}*/ 00186 }; 00187 00188 00216 class BCAdmin { 00217 public: 00218 typedef BCondition * bc_ptr_t; 00219 00220 private: 00221 std::map<int, bc_ptr_t, std::less<int> > _map; 00222 00223 public: 00228 template < 00229 class SP, 00230 class MAT, 00231 class XVEC, 00232 class BVEC 00233 > 00234 void apply(const SP& sp, 00235 MAT& A, 00236 XVEC& u, 00237 BVEC& f, 00238 bool preserve_symmetry = true) 00239 { 00240 u_int n_dof = sp.n_dof(); 00241 const SparsityPattern& spA = A.get_sparsity_pattern(); 00242 const std::size_t * rowstart = spA.get_rowstart_indices(); 00243 const u_int * colnum = spA.get_column_numbers(); 00244 for (u_int i = 0;i < n_dof;++ i) { 00245 int bm = sp.dofInfo(i).boundary_mark; 00246 if (bm == 0) continue; 00247 const bc_ptr_t bc = this->find(bm); 00249 if (bc == NULL) continue; 00250 if (bc->type() != BCondition::DIRICHLET) continue; 00251 bc->value((const void *)(&sp.dofInfo(i).interp_point), 00252 (void *)(&u(i))); 00253 00254 f(i) = A.diag_element(i)*u(i); 00255 for (u_int j = rowstart[i] + 1;j < rowstart[i + 1];++ j) { 00256 A.global_entry(j) -= A.global_entry(j); 00257 } 00258 if (preserve_symmetry) { 00259 for (u_int j = rowstart[i] + 1;j < rowstart[i + 1];++ j) { 00260 u_int k = colnum[j]; 00261 const u_int * p = std::find(&colnum[rowstart[k] + 1], 00262 &colnum[rowstart[k + 1]], i); 00263 if (p != &colnum[rowstart[k+1]]) { 00264 u_int l = p - &colnum[rowstart[0]]; 00265 f(k) -= A.global_entry(l)*f(i)/A.diag_element(i); 00266 A.global_entry(l) -= A.global_entry(l); 00267 } 00268 } 00269 } 00270 } 00271 } 00272 00277 template < 00278 class SP, 00279 class VEC 00280 > 00281 void clear_entry(const SP& sp, 00282 VEC& f) 00283 { 00284 u_int n_dof = sp.n_dof(); 00285 for (u_int i = 0;i < n_dof;++ i) { 00286 int bm = sp.dof_info(i).bound_mark(); 00287 if (bm == 0) continue; 00288 const bc_ptr_t bc = find(bm); 00289 if (bc != NULL) f(i) -= f(i); 00290 } 00291 } 00292 00299 template <class BC> 00300 void add(BC& b) 00301 { 00302 const std::set<int, std::less<int> >& bm = b.bound_mark(); 00303 std::set<int, std::less<int> >::const_iterator 00304 the_bm = bm.begin(), end_bm = bm.end(); 00305 for (;the_bm != end_bm;++ the_bm) { 00306 _map[*the_bm] = &b; 00307 } 00308 } 00309 00315 bc_ptr_t find(int bm) const 00316 { 00317 std::map<int, bc_ptr_t, std::less<int> >::const_iterator 00318 the_ptr = _map.find(bm); 00319 if (the_ptr != _map.end()) { 00320 return the_ptr->second; 00321 } else { 00322 return NULL; 00323 } 00324 } 00325 }; 00326 00327 #endif // __BoundaryCondition_h__ 00328