/usr/lib/petscdir/3.1/include/petscmesh_solvers.hh is in libpetsc3.1-dev 3.1.dfsg-11ubuntu1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 | #if !defined(__PETSCMESH_SOLVERS_HH)
#define __PETSCMESH_SOLVERS_HH
#include <petscmesh.hh>
#include <petscpc.h>
using ALE::Obj;
template<typename Section, typename Order>
void constructFieldSplit(const Obj<Section>& section, const Obj<Order>& globalOrder, Vec v, PC fieldSplit) {
const typename Section::chart_type& chart = section->getChart();
PetscInt space = 0;
PetscErrorCode ierr;
PetscInt total = 0;
for(typename std::vector<Obj<typename Section::atlas_type> >::const_iterator s_iter = section->getSpaces().begin(); s_iter != section->getSpaces().end(); ++s_iter, ++space) {
PetscInt n = section->size(space);
std::cout << "Space " << space << ": size " << n << std::endl;
total += n;
}
PetscInt localSize;
VecGetLocalSize(v, &localSize);
std::cout << "Vector local size " << localSize << std::endl;
assert(localSize == total);
space = 0;
for(typename std::vector<Obj<typename Section::atlas_type> >::const_iterator s_iter = section->getSpaces().begin(); s_iter != section->getSpaces().end(); ++s_iter, ++space) {
PetscInt n = section->size(space);
PetscInt i = -1;
PetscInt *idx;
IS is;
ierr = PetscMalloc(n * sizeof(PetscInt), &idx);CHKERRXX(ierr);
for(typename Section::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
const int dim = section->getFiberDimension(*c_iter, space);
const int cDim = section->getConstraintDimension(*c_iter, space);
if (dim > cDim) {
int off = globalOrder->getIndex(*c_iter);
for(int s = 0; s < space; ++s) {
off += section->getConstrainedFiberDimension(*c_iter, s);
}
if (cDim) {
// This is potentially dangerous
// These constraints dofs are for SINGLE FIELDS, not the entire point (confusing)
const int *cDofs = section->getConstraintDof(*c_iter, space);
for(int d = 0, c = 0, k = 0; d < dim; ++d) {
if ((c < cDim) && (cDofs[c] == d)) {
std::cout << " Ignored " << (off+k) << " at local pos " << d << " for point " << (*c_iter) << std::endl;
++c;
continue;
}
idx[++i] = off+k;
std::cout << "Added " << (off+k) << " at pos " << i << " for point " << (*c_iter) << std::endl;
++k;
}
} else {
for(int d = 0; d < dim; ++d) {
idx[++i] = off+d;
std::cout << "Added " << (off+d) << " at pos " << i << " for point " << (*c_iter) << std::endl;
}
}
}
}
if (i != n-1) {throw PETSc::Exception("Invalid fibration numbering");}
ierr = ISCreateGeneralNC(section->comm(), n, idx, &is);CHKERRXX(ierr);
ierr = PCFieldSplitSetIS(fieldSplit, is);CHKERRXX(ierr);
}
};
#endif // __PETSCMESH_SOLVERS_HH
|