Line data Source code
1 : /*
2 : * Copyright (c) 2010-2015: G-CSC, Goethe University Frankfurt
3 : * Author: Andreas Vogel
4 : *
5 : * This file is part of UG4.
6 : *
7 : * UG4 is free software: you can redistribute it and/or modify it under the
8 : * terms of the GNU Lesser General Public License version 3 (as published by the
9 : * Free Software Foundation) with the following additional attribution
10 : * requirements (according to LGPL/GPL v3 §7):
11 : *
12 : * (1) The following notice must be displayed in the Appropriate Legal Notices
13 : * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 : *
15 : * (2) The following notice must be displayed at a prominent place in the
16 : * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 : *
18 : * (3) The following bibliography is recommended for citation and must be
19 : * preserved in all covered files:
20 : * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 : * parallel geometric multigrid solver on hierarchically distributed grids.
22 : * Computing and visualization in science 16, 4 (2013), 151-164"
23 : * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 : * flexible software system for simulating pde based models on high performance
25 : * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 : *
27 : * This program is distributed in the hope that it will be useful,
28 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 : * GNU Lesser General Public License for more details.
31 : */
32 :
33 : #ifndef __H__UG__LIB_DISC__COMMON__LOCAL_ALGEBRA__
34 : #define __H__UG__LIB_DISC__COMMON__LOCAL_ALGEBRA__
35 :
36 : //#define UG_LOCALALGEBRA_ASSERT(cond, exp)
37 : // include define below to assert arrays used in stabilization
38 : #define UG_LOCALALGEBRA_ASSERT(cond, exp) UG_ASSERT((cond), exp)
39 :
40 : #include <vector>
41 :
42 : #include "./multi_index.h"
43 : #include "./function_group.h"
44 : #include "lib_algebra/small_algebra/small_algebra.h"
45 :
46 : namespace ug{
47 :
48 :
49 0 : class LocalIndices
50 : {
51 : public:
52 : /// Index type used by algebra
53 : typedef size_t index_type;
54 :
55 : /// Component type used by algebra
56 : typedef index_type comp_type;
57 :
58 : public:
59 : /// Default Constructor
60 : LocalIndices() {};
61 :
62 : /// sets the number of functions
63 : void resize_fct(size_t numFct)
64 : {
65 0 : m_vvIndex.resize(numFct);
66 0 : m_vLFEID.resize(numFct);
67 0 : }
68 :
69 : /// sets the local finite element id for a function
70 : void set_lfeID(size_t fct, const LFEID& lfeID)
71 : {
72 : UG_ASSERT(fct < m_vLFEID.size(), "Invalid index: "<<fct);
73 : m_vLFEID[fct] = lfeID;
74 : }
75 :
76 : /// returns the local finite element id of a function
77 : const LFEID& local_finite_element_id(size_t fct) const
78 : {
79 : UG_ASSERT(fct < m_vLFEID.size(), "Invalid index: "<<fct);
80 : return m_vLFEID[fct];
81 : }
82 :
83 : /// sets the number of dofs of a function
84 : void resize_dof(size_t fct, size_t numDoF)
85 : {
86 : check_fct(fct);
87 0 : m_vvIndex[fct].resize(numDoF);
88 : }
89 :
90 : /// clears the dofs of a function
91 : void clear_dof(size_t fct) {resize_dof(fct, 0);}
92 :
93 : /// reserves memory for the number of dofs
94 : void reserve_dof(size_t fct, size_t numDoF)
95 : {
96 : check_fct(fct);
97 : m_vvIndex[fct].reserve(numDoF);
98 : }
99 :
100 : /// adds an index (increases size)
101 0 : void push_back_index(size_t fct, size_t index) {push_back_multi_index(fct,index,0);}
102 :
103 : /// adds an index (increases size)
104 : void push_back_multi_index(size_t fct, size_t index, size_t comp)
105 : {
106 : check_fct(fct);
107 0 : m_vvIndex[fct].push_back(DoFIndex(index,comp));
108 0 : }
109 :
110 : /// clears all fct
111 : void clear() {m_vvIndex.clear();}
112 :
113 : /// number of functions
114 : size_t num_fct() const {return m_vvIndex.size();}
115 :
116 : /// number of dofs for accessible function
117 : size_t num_dof(size_t fct) const
118 : {
119 : check_fct(fct);
120 : return m_vvIndex[fct].size();
121 : }
122 :
123 : /// number of dofs of all accessible (sum)
124 : size_t num_dof() const
125 : {
126 : size_t num = 0;
127 : for(size_t fct = 0; fct < num_fct(); ++fct) num += num_dof(fct);
128 : return num;
129 : }
130 :
131 : /// global algebra multi-index for (fct, dof)
132 : const DoFIndex& multi_index(size_t fct, size_t dof) const
133 : {
134 : check_dof(fct, dof);
135 : return m_vvIndex[fct][dof];
136 : }
137 :
138 : /// global algebra index for (fct, dof)
139 : index_type index(size_t fct, size_t dof) const
140 : {
141 : check_dof(fct, dof);
142 0 : return m_vvIndex[fct][dof][0];
143 : }
144 :
145 : /// global algebra index for (fct, dof)
146 : index_type& index(size_t fct, size_t dof)
147 : {
148 : check_dof(fct, dof);
149 : return m_vvIndex[fct][dof][0];
150 : }
151 :
152 : /// algebra comp for (fct, dof)
153 : comp_type comp(size_t fct, size_t dof) const
154 : {
155 : check_dof(fct, dof);
156 0 : return m_vvIndex[fct][dof][1];
157 : }
158 :
159 : /// algebra comp for (fct, dof)
160 : comp_type& comp(size_t fct, size_t dof)
161 : {
162 : check_dof(fct, dof);
163 : return m_vvIndex[fct][dof][1];
164 : }
165 :
166 : /// checks if the local index object references a given index
167 : bool contains_index(index_type ind)
168 : {
169 : for(size_t fct = 0; fct < num_fct(); fct++)
170 : for(size_t dof = 0; dof < num_dof(fct); dof++)
171 : if(m_vvIndex[fct][dof][0] == ind)
172 : return true;
173 : return false;
174 : }
175 :
176 : protected:
177 : /// checks correct fct index in debug mode
178 : inline void check_fct(size_t fct) const
179 : {
180 : UG_LOCALALGEBRA_ASSERT(fct < num_fct(), "Wrong index.");
181 : }
182 : /// checks correct dof index in debug mode
183 : inline void check_dof(size_t fct, size_t dof) const
184 : {
185 : check_fct(fct);
186 : UG_LOCALALGEBRA_ASSERT(dof < num_dof(fct), "Wrong index.");
187 : }
188 :
189 : protected:
190 : // Mapping (fct, dof) -> local index
191 : std::vector<std::vector<DoFIndex> > m_vvIndex;
192 :
193 : // Local finite element ids
194 : std::vector<LFEID> m_vLFEID;
195 : };
196 :
197 0 : class LocalVector
198 : {
199 : public:
200 : /// own type
201 : typedef LocalVector this_type;
202 :
203 : /// Type to store DoF values
204 : typedef number value_type;
205 :
206 : public:
207 : /// default Constructor
208 0 : LocalVector() : m_pIndex(NULL) {m_vvValue.clear();}
209 :
210 : /// Constructor
211 : LocalVector(const LocalIndices& ind) {resize(ind);}
212 :
213 : /// resize for current local indices
214 0 : void resize(const LocalIndices& ind)
215 : {
216 0 : m_pIndex = &ind;
217 0 : m_vvValue.resize(ind.num_fct());
218 0 : m_vvValueAcc.resize(m_vvValue.size());
219 0 : for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
220 0 : m_vvValue[fct].resize(ind.num_dof(fct));
221 0 : access_all();
222 0 : }
223 :
224 : /// get current local indices
225 0 : const LocalIndices& get_indices() const {return *m_pIndex;}
226 :
227 : ///////////////////////////
228 : // vector functions
229 : ///////////////////////////
230 :
231 0 : this_type& operator=(const this_type& other)
232 : {
233 0 : m_pIndex = other.m_pIndex;
234 : const size_t numFcts = m_pIndex->num_fct();
235 0 : m_vvValue.resize(numFcts);
236 0 : for (size_t fct = 0; fct < numFcts; ++fct)
237 0 : m_vvValue[fct].resize(m_pIndex->num_dof(fct));
238 0 : m_vvValueAcc.resize(numFcts);
239 0 : if (other.m_pFuncMap)
240 : access_by_map(*other.m_pFuncMap);
241 : else
242 0 : access_all();
243 :
244 0 : return *this;
245 : }
246 :
247 : /// set all components of the vector
248 0 : this_type& operator=(number val)
249 : {
250 0 : for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
251 0 : for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
252 0 : m_vvValue[fct][dof] = val;
253 0 : return *this;
254 : }
255 :
256 : /// multiply all components of the vector
257 : this_type& operator*(number val)
258 : {
259 : return this->operator*=(val);
260 : }
261 :
262 : /// multiply all components of the vector
263 : this_type& operator*=(number val)
264 : {
265 : for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
266 : for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
267 : m_vvValue[fct][dof] *= val;
268 : return *this;
269 : }
270 :
271 : /// add a local vector
272 : this_type& operator+=(const this_type& rhs)
273 : {
274 : UG_LOCALALGEBRA_ASSERT(m_pIndex==rhs.m_pIndex, "Not same indices.");
275 : for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
276 : for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
277 : m_vvValue[fct][dof] += rhs.m_vvValue[fct][dof];
278 : return *this;
279 : }
280 :
281 : /// subtract a local vector
282 : this_type& operator-=(const this_type& rhs)
283 : {
284 : UG_LOCALALGEBRA_ASSERT(m_pIndex==rhs.m_pIndex, "Not same indices.");
285 : for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
286 : for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
287 : m_vvValue[fct][dof] -= rhs.m_vvValue[fct][dof];
288 : return *this;
289 : }
290 :
291 : /// add a scaled vector
292 0 : this_type& scale_append(number s, const this_type& rhs)
293 : {
294 : UG_LOCALALGEBRA_ASSERT(m_pIndex==rhs.m_pIndex, "Not same indices.");
295 0 : for(size_t fct = 0; fct < m_vvValue.size(); ++fct)
296 0 : for(size_t dof = 0; dof < m_vvValue[fct].size(); ++dof)
297 0 : m_vvValue[fct][dof] += s * rhs.m_vvValue[fct][dof];
298 0 : return *this;
299 : }
300 :
301 : ///////////////////////////
302 : // restricted DoF access
303 : ///////////////////////////
304 :
305 : /// access only part of the functions using mapping (restrict functions)
306 : void access_by_map(const FunctionIndexMapping& funcMap)
307 : {
308 0 : m_pFuncMap = &funcMap;
309 0 : for(size_t i = 0; i < funcMap.num_fct(); ++i)
310 : {
311 : const size_t mapFct = funcMap[i];
312 0 : m_vvValueAcc[i] = &(m_vvValue[mapFct][0]);
313 : }
314 : }
315 :
316 : /// access all functions
317 0 : void access_all()
318 : {
319 0 : m_pFuncMap = NULL;
320 :
321 0 : if(m_pIndex==NULL) {m_vvValueAcc.clear(); return;}
322 :
323 0 : for(size_t i = 0; i < m_pIndex->num_fct(); ++i)
324 0 : m_vvValueAcc[i] = &(m_vvValue[i][0]);
325 : }
326 :
327 : /// returns the number of currently accessible functions
328 : size_t num_fct() const
329 : {
330 : if(m_pFuncMap == NULL) return m_vvValue.size();
331 : return m_pFuncMap->num_fct();
332 : }
333 :
334 : /// returns the local finite element id of a function
335 : const LFEID& local_finite_element_id(size_t fct) const
336 : {
337 : UG_ASSERT(m_pIndex != NULL, "No indices present");
338 : check_fct(fct);
339 : if(m_pFuncMap == NULL) return m_pIndex->local_finite_element_id(fct);
340 : else return m_pIndex->local_finite_element_id((*m_pFuncMap)[fct]);
341 : }
342 :
343 : /// returns the number of dofs for the currently accessible function
344 : size_t num_dof(size_t fct) const
345 : {
346 : check_fct(fct);
347 : if(m_pFuncMap == NULL) return m_vvValue[fct].size();
348 : else return m_vvValue[ (*m_pFuncMap)[fct] ].size();
349 : }
350 :
351 : /// access to dof of currently accessible function fct
352 : number& operator()(size_t fct, size_t dof)
353 : {
354 : check_dof(fct,dof);
355 0 : return m_vvValueAcc[fct][dof];
356 : }
357 :
358 : /// const access to dof of currently accessible function fct
359 : number operator()(size_t fct, size_t dof) const
360 : {
361 : check_dof(fct,dof);
362 0 : return m_vvValueAcc[fct][dof];
363 : }
364 :
365 : ///////////////////////////
366 : // all DoF access
367 : ///////////////////////////
368 :
369 : /// returns the number of all functions
370 : size_t num_all_fct() const {return m_vvValue.size();}
371 :
372 : /// returns the number of dofs for a function (unrestricted functions)
373 : size_t num_all_dof(size_t fct) const {check_all_fct(fct); return m_vvValue[fct].size();}
374 :
375 : /// access to dof of a fct (unrestricted functions)
376 : number& value(size_t fct, size_t dof){check_all_dof(fct,dof);return m_vvValue[fct][dof];}
377 :
378 : /// const access to dof of a fct (unrestricted functions)
379 : const number& value(size_t fct, size_t dof) const{check_all_dof(fct,dof);return m_vvValue[fct][dof];}
380 :
381 : protected:
382 : /// checks correct fct index in debug mode
383 : inline void check_fct(size_t fct) const
384 : {
385 : UG_LOCALALGEBRA_ASSERT(fct < num_fct(), "Wrong index: fct: "<<fct<<
386 : " of num_fct: "<<num_fct());
387 : }
388 : /// checks correct dof index in debug mode
389 : inline void check_dof(size_t fct, size_t dof) const
390 : {
391 : check_fct(fct);
392 : UG_LOCALALGEBRA_ASSERT(dof < num_dof(fct), "Wrong index: dof: "
393 : <<dof<<", num_dof: "<<num_dof(fct)<<" of fct: "<<fct);
394 : }
395 : /// checks correct fct index in debug mode
396 : inline void check_all_fct(size_t fct) const
397 : {
398 : UG_LOCALALGEBRA_ASSERT(fct < num_all_fct(), "Wrong index.");
399 : }
400 : /// checks correct dof index in debug mode
401 : inline void check_all_dof(size_t fct, size_t dof) const
402 : {
403 : check_all_fct(fct);
404 : UG_LOCALALGEBRA_ASSERT(dof < num_all_dof(fct), "Wrong index.");
405 : }
406 :
407 : protected:
408 : /// Indices
409 : const LocalIndices* m_pIndex;
410 :
411 : /// Access Mapping
412 : const FunctionIndexMapping* m_pFuncMap;
413 :
414 : /// Entries (fct, dof)
415 : std::vector<value_type*> m_vvValueAcc;
416 :
417 : /// Entries (fct, dof)
418 : std::vector<std::vector<value_type> > m_vvValue;
419 : };
420 :
421 : class LocalMatrix
422 : {
423 : public:
424 : /// own type
425 : typedef LocalMatrix this_type;
426 :
427 : /// Entry type used by algebra
428 : typedef number value_type;
429 :
430 : public:
431 : /// Constructor
432 : LocalMatrix() :
433 0 : m_pRowIndex(NULL), m_pColIndex(NULL) ,
434 0 : m_pRowFuncMap(NULL), m_pColFuncMap(NULL)
435 : {}
436 :
437 : /// Constructor
438 : LocalMatrix(const LocalIndices& rowInd, const LocalIndices& colInd)
439 : : m_pRowFuncMap(NULL), m_pColFuncMap(NULL)
440 : {
441 : resize(rowInd, colInd);
442 : }
443 :
444 : /// resize for same ind in row and column
445 0 : void resize(const LocalIndices& ind) {resize(ind, ind);}
446 :
447 : /// resize for current local indices
448 0 : void resize(const LocalIndices& rowInd, const LocalIndices& colInd)
449 : {
450 0 : m_pRowIndex = &rowInd;
451 0 : m_pColIndex = &colInd;
452 :
453 0 : m_CplMat.resize(rowInd.num_fct(), colInd.num_fct());
454 0 : m_CplMatAcc.resize(rowInd.num_fct(), colInd.num_fct());
455 :
456 0 : for(size_t fct1 = 0; fct1 < m_CplMat.num_rows(); ++fct1)
457 0 : for(size_t fct2 = 0; fct2 < m_CplMat.num_cols(); ++fct2)
458 : {
459 0 : m_CplMat(fct1, fct2).resize(colInd.num_dof(fct1),
460 : rowInd.num_dof(fct2));
461 : }
462 :
463 0 : access_all();
464 0 : }
465 :
466 : /// get current local indices
467 0 : const LocalIndices& get_row_indices() const {return *m_pRowIndex;}
468 :
469 : /// get current local indices
470 0 : const LocalIndices& get_col_indices() const {return *m_pColIndex;}
471 :
472 : ///////////////////////////
473 : // Matrix functions
474 : ///////////////////////////
475 :
476 : /// set all entries
477 0 : this_type& operator=(number val)
478 : {
479 0 : for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
480 0 : for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
481 0 : m_CplMat(fct1,fct2) = val;
482 0 : return *this;
483 : }
484 :
485 : /// multiply all entries
486 : this_type& operator*(number val)
487 : {
488 : return this->operator*=(val);
489 : }
490 :
491 : /// multiply matrix
492 0 : this_type& operator*=(number val)
493 : {
494 0 : for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
495 0 : for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
496 : m_CplMat(fct1,fct2) *= val;
497 0 : return *this;
498 : }
499 :
500 : /// add matrix
501 : this_type& operator+=(const this_type& rhs)
502 : {
503 : UG_LOCALALGEBRA_ASSERT(m_pRowIndex==rhs.m_pRowIndex &&
504 : m_pColIndex==rhs.m_pColIndex, "Not same indices.");
505 : for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
506 : for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
507 : m_CplMat(fct1,fct2) += rhs.m_CplMat(fct1,fct2);
508 : return *this;
509 : }
510 :
511 : /// subtract matrix
512 : this_type& operator-=(const this_type& rhs)
513 : {
514 : UG_LOCALALGEBRA_ASSERT(m_pRowIndex==rhs.m_pRowIndex &&
515 : m_pColIndex==rhs.m_pColIndex, "Not same indices.");
516 : for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
517 : for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
518 : m_CplMat(fct1,fct2) -= rhs.m_CplMat(fct1,fct2);
519 : return *this;
520 : }
521 :
522 : /// add scaled matrix
523 0 : this_type& scale_append(number s, const this_type& rhs)
524 : {
525 : UG_LOCALALGEBRA_ASSERT(m_pRowIndex==rhs.m_pRowIndex &&
526 : m_pColIndex==rhs.m_pColIndex, "Not same indices.");
527 0 : for(size_t fct1=0; fct1 < m_CplMat.num_rows(); ++fct1)
528 0 : for(size_t fct2=0; fct2 < m_CplMat.num_cols(); ++fct2)
529 0 : MatScaleAppend(m_CplMat(fct1,fct2), s, rhs.m_CplMat(fct1,fct2));
530 0 : return *this;
531 : }
532 :
533 : ///////////////////////////
534 : // restricted DoF access
535 : ///////////////////////////
536 :
537 : /// access only part of the functions using mapping (restrict functions)
538 : void access_by_map(const FunctionIndexMapping& funcMap)
539 : {
540 0 : access_by_map(funcMap, funcMap);
541 : }
542 :
543 : /// access only part of the functions using mapping (restrict functions)
544 0 : void access_by_map(const FunctionIndexMapping& rowFuncMap,
545 : const FunctionIndexMapping& colFuncMap)
546 : {
547 0 : m_pRowFuncMap = &rowFuncMap;
548 0 : m_pColFuncMap = &colFuncMap;
549 :
550 0 : for(size_t i = 0; i < m_pRowFuncMap->num_fct(); ++i)
551 0 : for(size_t j = 0; j < m_pColFuncMap->num_fct(); ++j)
552 : {
553 : const size_t rowMapFct = rowFuncMap[i];
554 : const size_t colMapFct = colFuncMap[j];
555 :
556 0 : m_CplMatAcc(i,j) = &(m_CplMat(rowMapFct, colMapFct));
557 : }
558 0 : }
559 :
560 : /// access all functions
561 0 : void access_all()
562 : {
563 0 : m_pRowFuncMap = NULL;
564 0 : m_pColFuncMap = NULL;
565 :
566 0 : if(m_pRowIndex==NULL) {m_CplMatAcc.resize(0,0); return;}
567 :
568 0 : for(size_t i = 0; i < m_pRowIndex->num_fct(); ++i)
569 0 : for(size_t j = 0; j < m_pColIndex->num_fct(); ++j)
570 0 : m_CplMatAcc(i,j) = &(m_CplMat(i,j));
571 : }
572 :
573 : /// returns the number of currently accessible (restricted) functions
574 : size_t num_row_fct() const
575 : {
576 : if(m_pRowFuncMap != NULL) return m_pRowFuncMap->num_fct();
577 : return m_CplMatAcc.num_rows();
578 : }
579 :
580 : /// returns the number of currently accessible (restricted) functions
581 : size_t num_col_fct() const
582 : {
583 : if(m_pColFuncMap != NULL) return m_pColFuncMap->num_fct();
584 : return m_CplMatAcc.num_cols();
585 : }
586 :
587 : /// returns the number of dofs for the currently accessible (restricted) function
588 : size_t num_row_dof(size_t fct) const
589 : {
590 : if(m_CplMat.num_rows()==0) return 0;
591 : if(m_pRowFuncMap == NULL) return m_CplMat(fct, 0).num_rows();
592 : else return m_CplMat( (*m_pRowFuncMap)[fct], 0).num_rows();
593 : }
594 :
595 : /// returns the number of dofs for the currently accessible (restricted) function
596 : size_t num_col_dof(size_t fct) const
597 : {
598 : if(m_CplMat.num_cols()==0) return 0;
599 : if(m_pRowFuncMap == NULL) return m_CplMat(0, fct).num_cols();
600 : else return m_CplMat( 0, (*m_pColFuncMap)[fct]).num_cols();
601 : }
602 :
603 : /// access to (restricted) coupling (rowFct, rowDoF) x (colFct, colDoF)
604 : number& operator()(size_t rowFct, size_t rowDoF,
605 : size_t colFct, size_t colDoF)
606 : {
607 : check_dof(rowFct, rowDoF, colFct, colDoF);
608 0 : return (*m_CplMatAcc(rowFct,colFct))(rowDoF, colDoF);
609 : }
610 :
611 : /// const access to (restricted) coupling (rowFct, rowDoF) x (colFct, colDoF)
612 : number operator()(size_t rowFct, size_t rowDoF,
613 : size_t colFct, size_t colDoF) const
614 : {
615 : check_dof(rowFct, rowDoF, colFct, colDoF);
616 : return (*m_CplMatAcc(rowFct,colFct))(rowDoF, colDoF);
617 : }
618 :
619 : ///////////////////////////
620 : // all DoF access
621 : ///////////////////////////
622 :
623 : /// returns the number of all functions
624 : size_t num_all_row_fct() const{ return m_CplMat.num_rows();}
625 :
626 : /// returns the number of all functions
627 : size_t num_all_col_fct() const{return m_CplMat.num_cols();}
628 :
629 : /// returns the number of dofs for a function
630 : size_t num_all_row_dof(size_t fct) const {return m_CplMat(fct, 0).num_rows();}
631 :
632 : /// returns the number of dofs for a function
633 : size_t num_all_col_dof(size_t fct) const {return m_CplMat(0, fct).num_cols();}
634 :
635 : /// access to coupling (rowFct, rowDoF) x (colFct, colDoF)
636 : number& value(size_t rowFct, size_t rowDoF,
637 : size_t colFct, size_t colDoF)
638 : {
639 : check_all_dof(rowFct, rowDoF, colFct, colDoF);
640 : return (m_CplMat(rowFct,colFct))(rowDoF, colDoF);
641 : }
642 :
643 : /// const access to coupling (rowFct, rowDoF) x (colFct, colDoF)
644 : number value(size_t rowFct, size_t rowDoF,
645 : size_t colFct, size_t colDoF) const
646 : {
647 : check_all_dof(rowFct, rowDoF, colFct, colDoF);
648 0 : return (m_CplMat(rowFct,colFct))(rowDoF, colDoF);
649 : }
650 :
651 : protected:
652 : /// checks correct (fct1,fct2) index in debug mode
653 : inline void check_fct(size_t rowFct, size_t colFct) const
654 : {
655 : UG_LOCALALGEBRA_ASSERT(rowFct < num_row_fct(), "Wrong index.");
656 : UG_LOCALALGEBRA_ASSERT(colFct < num_col_fct(), "Wrong index.");
657 : }
658 : /// checks correct (dof1,dof2) index in debug mode
659 : inline void check_dof(size_t rowFct, size_t rowDoF,
660 : size_t colFct, size_t colDoF) const
661 : {
662 : check_fct(rowFct, colFct);
663 : UG_LOCALALGEBRA_ASSERT(rowDoF < num_row_dof(rowFct), "Wrong index.");
664 : UG_LOCALALGEBRA_ASSERT(colDoF < num_col_dof(colFct), "Wrong index.");
665 : }
666 : /// checks correct (fct1,fct2) index in debug mode
667 : inline void check_all_fct(size_t rowFct, size_t colFct) const
668 : {
669 : UG_LOCALALGEBRA_ASSERT(rowFct < num_all_row_fct(), "Wrong index.");
670 : UG_LOCALALGEBRA_ASSERT(colFct < num_all_col_fct(), "Wrong index.");
671 : }
672 : /// checks correct (dof1,dof2) index in debug mode
673 : inline void check_all_dof(size_t rowFct, size_t rowDoF,
674 : size_t colFct, size_t colDoF) const
675 : {
676 : check_all_fct(rowFct, colFct);
677 : UG_LOCALALGEBRA_ASSERT(rowDoF < num_all_row_dof(rowFct), "Wrong index.");
678 : UG_LOCALALGEBRA_ASSERT(colDoF < num_all_col_dof(colFct), "Wrong index.");
679 : }
680 :
681 : protected:
682 : // Row indices
683 : const LocalIndices* m_pRowIndex;
684 :
685 : // Column indices
686 : const LocalIndices* m_pColIndex;
687 :
688 : /// Row Access Mapping
689 : const FunctionIndexMapping* m_pRowFuncMap;
690 :
691 : /// Column Access Mapping
692 : const FunctionIndexMapping* m_pColFuncMap;
693 :
694 : // \todo: Think of a better (faster) storage (one plain array?)
695 :
696 : // type of cpl matrices between two functions
697 : typedef DenseMatrix<VariableArray2<value_type> > LocalCplMatrix;
698 :
699 : // type of Func-Coupling matrices
700 : typedef DenseMatrix<VariableArray2<LocalCplMatrix> > FctCplMatrix;
701 :
702 : // type of Func-Coupling pointer matrices
703 : typedef DenseMatrix<VariableArray2<LocalCplMatrix*> > FctCplAccMatrix;
704 :
705 : // Entries (fct1, fct2, dof1, dof2)
706 : FctCplMatrix m_CplMat;
707 :
708 : // Entries (fct1, fct2, dof1, dof2)
709 : FctCplAccMatrix m_CplMatAcc;
710 : };
711 :
712 : inline
713 : std::ostream& operator<< (std::ostream& outStream, const ug::LocalMatrix& mat)
714 : {
715 : for(size_t fct1 = 0; fct1 < mat.num_row_fct(); ++fct1)
716 : for(size_t fct2 = 0; fct2 < mat.num_col_fct(); ++fct2)
717 : for(size_t dof1 = 0; dof1 < mat.num_row_dof(fct1); ++dof1)
718 : for(size_t dof2 = 0; dof2 < mat.num_col_dof(fct2); ++dof2)
719 : {
720 : outStream << "[("<< fct1 << "," << dof1 << ") x ("
721 : << fct2 << "," << dof2 << ")]: "
722 : << mat(fct1, dof1, fct2, dof2) << "\n";
723 : }
724 : return outStream;
725 : }
726 :
727 : inline
728 : std::ostream& operator<< (std::ostream& outStream, const ug::LocalVector& vec)
729 : {
730 : for(size_t fct = 0; fct < vec.num_fct(); ++fct)
731 : for(size_t dof = 0; dof < vec.num_dof(fct); ++dof)
732 : {
733 : outStream << "["<<fct<<","<<dof<<"]:" << vec(fct, dof) << "\n";
734 : }
735 : return outStream;
736 : }
737 :
738 : template <typename TVector>
739 0 : void GetLocalVector(LocalVector& lvec, const TVector& vec)
740 : {
741 : const LocalIndices& ind = lvec.get_indices();
742 :
743 0 : for(size_t fct=0; fct < lvec.num_all_fct(); ++fct)
744 0 : for(size_t dof=0; dof < lvec.num_all_dof(fct); ++dof)
745 : {
746 : const size_t index = ind.index(fct,dof);
747 : const size_t comp = ind.comp(fct,dof);
748 0 : lvec.value(fct,dof) = BlockRef(vec[index], comp);
749 : }
750 0 : }
751 :
752 : template <typename TVector>
753 0 : void AddLocalVector(TVector& vec, const LocalVector& lvec)
754 : {
755 : const LocalIndices& ind = lvec.get_indices();
756 :
757 0 : for(size_t fct=0; fct < lvec.num_all_fct(); ++fct)
758 0 : for(size_t dof=0; dof < lvec.num_all_dof(fct); ++dof)
759 : {
760 : const size_t index = ind.index(fct,dof);
761 : const size_t comp = ind.comp(fct,dof);
762 0 : BlockRef(vec[index], comp) += lvec.value(fct,dof);
763 : }
764 0 : }
765 :
766 : template <typename TMatrix>
767 0 : void AddLocalMatrixToGlobal(TMatrix& mat, const LocalMatrix& lmat)
768 : {
769 : //PROFILE_FUNC_GROUP("algebra")
770 : const LocalIndices& rowInd = lmat.get_row_indices();
771 : const LocalIndices& colInd = lmat.get_col_indices();
772 :
773 0 : for(size_t fct1=0; fct1 < lmat.num_all_row_fct(); ++fct1)
774 0 : for(size_t dof1=0; dof1 < lmat.num_all_row_dof(fct1); ++dof1)
775 : {
776 : const size_t rowIndex = rowInd.index(fct1,dof1);
777 : const size_t rowComp = rowInd.comp(fct1,dof1);
778 :
779 0 : for(size_t fct2=0; fct2 < lmat.num_all_col_fct(); ++fct2)
780 0 : for(size_t dof2=0; dof2 < lmat.num_all_col_dof(fct2); ++dof2)
781 : {
782 : const size_t colIndex = colInd.index(fct2,dof2);
783 : const size_t colComp = colInd.comp(fct2,dof2);
784 :
785 : BlockRef(mat(rowIndex, colIndex), rowComp, colComp)
786 0 : += lmat.value(fct1,dof1,fct2,dof2);
787 : }
788 : }
789 0 : }
790 :
791 : } // end namespace ug
792 :
793 : #endif /* __H__UG__LIB_DISC__COMMON__LOCAL_ALGEBRA__ */
|
