3 #ifndef DUNE_PETSCNESTEDVECTORBACKEND_HH
4 #define DUNE_PETSCNESTEDVECTORBACKEND_HH
10 #include<dune/common/fvector.hh>
22 class PetscNestedVectorBackend;
24 class PetscNestedVectorContainer
27 typedef std::vector<shared_ptr<PetscVectorContainer> > ChildContainer;
30 typedef PetscScalar ElementType;
31 typedef ElementType E;
32 typedef Vec ContainerType;
35 typedef block_type field_type;
37 typedef std::size_t size_type;
38 typedef PetscVectorBackend Backend;
40 struct create_child_vectors
41 :
public TypeTree::DirectChildrenVisitor
42 ,
public TypeTree::DynamicTraversal
44 create_child_vectors(ChildContainer& children)
48 template<
typename T,
typename Child,
typename TreePath,
typename ChildIndex>
49 void beforeChild(T&& t, Child&& child, TreePath treePath, ChildIndex childIndex)
51 _children[childIndex] = make_shared<PetscVectorContainer>(child);
54 ChildContainer& _children;
57 template<
typename GFS>
58 PetscNestedVectorContainer (
const GFS& gfs)
59 : _data(gfs.globalSize())
60 , _sub_data(GFS::CHILDREN)
63 const size_t _blockCount = GFS::CHILDREN;
64 ChildContainer children(_blockCount);
65 create_child_vectors ccv(children);
66 TypeTree::applyToTree(gfs,ccv);
67 Vec child_array[_blockCount];
68 for (
size_t i = 0; i < _blockCount; ++i)
69 child_array[i] = children[i]->_v;
70 PETSC_CALL(VecCreateNest(PETSC_COMM_SELF,_blockCount,PETSC_NULL,child_array,&_v));
73 template<
typename GFS>
74 PetscNestedVectorContainer (
const GFS& gfs,
const E&
e)
75 : _data(gfs.globalSize())
76 , _sub_data(GFS::CHILDREN)
79 const size_t _blockCount = GFS::CHILDREN;
80 ChildContainer children(_blockCount);
81 create_child_vectors ccv(children);
82 TypeTree::applyToTree(gfs,ccv);
83 Vec child_array[_blockCount];
84 for (
size_t i = 0; i < _blockCount; ++i)
85 child_array[i] = *(children[i]);
86 PETSC_CALL(VecCreateNest(PETSC_COMM_SELF,_blockCount,PETSC_NULL,child_array,&_v));
90 PetscNestedVectorContainer (
const PetscNestedVectorContainer& rhs)
91 : _data(rhs._data.size())
92 , _sub_data(rhs._sub_data.size())
96 PETSC_CALL(VecDuplicate(rhs._v,&_v));
97 PETSC_CALL(VecCopy(rhs._v,_v));
100 PetscNestedVectorContainer& operator= (
const PetscNestedVectorContainer& rhs)
104 PETSC_CALL(VecCopy(rhs._v,_v));
108 ~PetscNestedVectorContainer()
111 PETSC_CALL(VecDestroy(&_v));
117 PETSC_CALL(VecGetLocalSize(_v,&size));
122 PetscNestedVectorContainer& operator= (
const E& e)
125 PETSC_CALL(VecSet(_v,e));
129 PetscNestedVectorContainer& operator*= (
const E& e)
132 PETSC_CALL(VecScale(_v,e));
137 PetscNestedVectorContainer& operator+= (
const E& e)
140 PETSC_CALL(VecShift(_v,e));
144 PetscNestedVectorContainer& operator+= (
const PetscNestedVectorContainer& rhs)
146 return axpy(1.0,rhs);
149 PetscNestedVectorContainer& operator-= (
const PetscNestedVectorContainer& rhs)
151 return axpy(-1.0,rhs);
154 field_type& operator[](std::size_t i)
160 const field_type& operator[](std::size_t i)
const
170 PETSC_CALL(VecNorm(_v,NORM_2,&norm));
178 PETSC_CALL(VecNorm(_v,NORM_1,&norm));
182 E infinity_norm()
const
186 PETSC_CALL(VecNorm(_v,NORM_INFINITY,&norm));
190 E operator*(
const PetscNestedVectorContainer& y)
const
193 PetscScalar dotproduct = 0;
194 PETSC_CALL(VecTDot(_v,y._v,&dotproduct));
198 E dot(
const PetscNestedVectorContainer& y)
const
201 PetscScalar dotproduct = 0;
202 PETSC_CALL(VecDot(_v,y._v,&dotproduct));
207 PetscNestedVectorContainer& axpy(
const E& a,
const PetscNestedVectorContainer& y)
211 PETSC_CALL(VecAXPY(_v,a,y._v));
216 ContainerType& base ()
222 const ContainerType& base ()
const
228 operator ContainerType&()
234 operator const ContainerType&()
const
264 size_t flatsize()
const
270 void std_copy_to (std::vector<X>& x)
const
273 std::copy(_data.begin(),_data.end(),x.begin());
277 void std_copy_from (
const std::vector<X>& x)
280 std::copy(x.begin(),x.end(),_data.begin());
283 PetscVectorContainer subVector(size_type i)
287 PETSC_CALL(VecNestGetSubVec(_v,i,&sub_vec));
288 return PetscVectorContainer(sub_vec,
false);
293 mutable std::vector<PetscScalar> _data;
294 mutable std::vector<PetscScalar*> _sub_data;
295 mutable bool _checkedIn;
303 PETSC_CALL(VecNestGetSubVecs(_v,&N,&sv));
304 auto it = _data.begin();
305 for (
int n = 0; n < N; ++n)
308 PETSC_CALL(VecGetLocalSize(sv[n],&CN));
309 auto endit = it + CN;
310 for(
auto childit = _sub_data[n]; it != endit; ++it, ++childit)
312 PETSC_CALL(VecRestoreArray(sv[n],&(_sub_data[n])));
318 void checkout()
const
324 PETSC_CALL(VecNestGetSubVecs(_v,&N,&sv));
325 auto it = _data.begin();
326 for (
int n = 0; n < N; ++n)
329 PETSC_CALL(VecGetArray(sv[n],&(_sub_data[n])));
330 PETSC_CALL(VecGetLocalSize(sv[n],&CN));
331 auto endit = it + CN;
332 for(
auto childit = _sub_data[n]; it != endit; ++it, ++childit)
344 class PetscNestedVectorBackend
361 typedef typename C::ElementType Type;
365 typedef typename std::size_t size_type;
370 static const typename C::field_type& access (
const C& c, size_type i)
378 static typename C::field_type& access (C& c, size_type i)
384 template<
typename T,
typename E>
385 struct BackendVectorSelectorHelper<PetscNestedVectorBackend,T,E>
387 typedef PetscNestedVectorContainer Type;
397 #endif // DUNE_PETSCNESTEDVECTORBACKEND_HH
const E & e
Definition: interpolate.hh:172