Thanks Timo,
This was very helpful, but now I am confused about how to go about
doing what I want.
So, I have a bunch of
typedef std::vector < PETScWrappers::MPI::BlockVector > Vector;
objects that are initialized via an fesystem that is indexed via
FEValuesExtractors into scalar and vector components. All I want to
be able to do with these Vectors is restrict to a .block() (say the
first scalar component of the fesystem) in order to perform some
matrix algebra on the subsystem (some scalar component) using petsc
solvers, etc. I couldn't see how to do this without turning them into
block vectors. So I tried to make them into block vectors:
typedef std::vector < PETScWrappers::MPI::BlockVector > SolutionVector;
But that doesn't work because now each block gets treated as the full
fesystem when I initialize using the dofhandler created using
something like:
for (unsigned int component=0; component< alphadim; ++component){
dof_handler.push_back (new DoFHandler<dim>(triangulation));
alpha.push_back (new FEValuesExtractors::Scalar(0));
sigma.push_back (new FEValuesExtractors::Vector(1));
fe_collection.push_back ( new FESystem<dim>(FE_DGQ<dim>(max_degree), 1,
FE_DGQ<dim>(max_degree),
dim));
quadrature_collection.push_back ( new QGauss<dim>(max_degree+2));
face_quadrature_collection.push_back ( new QGauss<dim-1>(max_degree+2));
}
So there are two dofhandlers in each block. Which is not what I want.
Ideally, all I want is to initialize the whole system as a petsc
vector (or blockvector) where the fesystem is associated to the blocks
appropriately, (like block(0) corresponds to the first scalar
extractor and block(1) corresponds to the first vector extractor, and
so on, in order). Surely there's an easy way to do this, but I am
having trouble figuring it out.
Many thanks again,
Evan
On Tue, Jul 3, 2012 at 9:04 AM, Timo Heister <[email protected]> wrote:
> Hi Evan,
>
>> ...
>> subdomain_solution[component].collect_sizes();
>> ...
>
> You should only call collect_sizes() at the end after you constructed
> all blocks of the BlockVector.
>
>> ...
>> void dealii::DataOut_DoFData<DH, patch_dim,
>> patch_space_dim>::add_data_vector(const VECTOR&, const
>> std::vector<std::string, std::allocator<std::string> >&,
>> ...
>> The violated condition was:
>> vec.size() == dofs->n_dofs()
>> Additional Information:
>> The vector has size 23064 but the DoFHandler objects says there are
>> 11532 degrees of freedom and there are 961 active cells.
>
> You create a BlockVector with two separate DoFHandlers but hand the
> whole BlockVector to DataOut::add_data_vector(). This function expects
> a vector size that fits to the one DoFHandler you handed to DataOut.
> You need to use my_vector.block(0) or something in add_data_vector().
>
>> But I am not sure how to do this while preserving the ghosting.
>> Unlike the trilinos blockvectors in step-32, the petsc blockvectors do
>> not seem to support indexset:
>
> They do it the way you did it above. I will try to add a similar
> function to PETSc:MPI:BlockVector that takes a std::vector<IndexSet>
> soon.
>
> Timo
>
> --
> Timo Heister
> http://www.math.tamu.edu/~heister/
_______________________________________________
dealii mailing list http://poisson.dealii.org/mailman/listinfo/dealii
