Re: [deal.II] error when trying to initialize PETSc preconditioners
Many thanks Timo. I tried PreconditionBlockJacobi and it worked. However, I am running a validation loop of progressively smaller meshes, and stops again with a error when I make the call to initialize for the second time which is expected since ''initialize()'' can only be called once. Of course, there a simple solution for this problem: declaring the preconditioner inside the loop rather than outside and pass it as an argument to the routine that solves the linear system. Yet, is there a better solution? Say for instance, as with trilinos preconditioners that have a ''clear()'' function member to releases the preconditioner from the matrix, releases memory and so on, so that you can reinitialize it with a new matrix (this is particularly useful in adaptivity loops, like in step-32). That's all for now. Many Thanks. Ignacio. On Friday, August 12, 2016 at 9:31:36 PM UTC-5, Timo Heister wrote: > > Hey Ignacio, > > good to hear from you! > > Sadly this is not documented well, but the ICC preconditioner is only > implemented for BAIJ matrices and not MPIAIJ (which we use even if you > run on one processor), see > https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCICC.html > > You should be able to make it work by using a > PETScWrappers::SparseMatrix instead of an > PETScWrappers::MPI::SparseMatrix. Or, probably a better idea, just use > a preconditioner that works with parallel matrices like > PreconditionBlockJacobi. > > Best, > Timo > > > On Fri, Aug 12, 2016 at 9:46 PM, Ignacio Tomas > wrote: > > Hi, everybody. I am just trying to initialize either ICC, Jacobi or > > ParaSalis preconditioner from PETSc, using a PETSc matrix in parallel. > So > > far I am running just one mpi process. After the call to ''.initialize'' > I > > get: > > > > > > > > > > [0]PETSC ERROR: No support for this operation for this object type > > [0]PETSC ERROR: Matrix format mpiaij does not have a built-in PETSc ICC > > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html > for > > trouble shooting. > > [0]PETSC ERROR: Petsc Release Version 3.5.3, Jan, 31, 2015 > > [0]PETSC ERROR: ./main on a arch-linux2-c-opt named nachote by nachotex > Fri > > Aug 12 20:28:49 2016 > > [0]PETSC ERROR: Configure options --download-fblaslapack=1 > > --with-shared-libraries=1 --download-hypre=1 --download-mumps=1 > > --download-spooles=1 --download-scalapack=1 --download-metis=1 > > --download-parmetis=1 --download-blacs=1 --with-debugging=0 --with-x=0 > > [0]PETSC ERROR: #1 MatGetFactor() line 3960 in > > > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/mat/interface/matrix.c > > [0]PETSC ERROR: #2 PCSetup_ICC() line 18 in > > > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/ksp/pc/impls/factor/icc/icc.c > > > > [0]PETSC ERROR: #3 PCSetUp() line 902 in > > > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/ksp/pc/interface/precon.c > > > terminate called after throwing an instance of > > 'dealii::LACExceptions::ExcPETScError' > > what(): > > > > An error occurred in line <384> of file > > > > > > > in function > > void dealii::PETScWrappers::PreconditionICC::initialize(const > > dealii::PETScWrappers::MatrixBase&, const > > dealii::PETScWrappers::PreconditionICC::AdditionalData&) > > The violated condition was: > > ierr == 0 > > The name and call sequence of the exception was: > > ExcPETScError(ierr) > > Additional Information: > > An error with error number 56 occurred while calling a PETSc function > > > > > > > > > > Any clue? Did I forget to install something? > > > > -- > > The deal.II project is located at http://www.dealii.org/ > > For mailing list/forum options, see > > https://groups.google.com/d/forum/dealii?hl=en > > --- > > You received this message because you are subscribed to the Google > Groups > > "deal.II User Group" group. > > To unsubscribe from this group and stop receiving emails from it, send > an > > email to dealii+un...@googlegroups.com . > > For more options, visit https://groups.google.com/d/optout. > > > > -- > Timo Heister > http://www.math.clemson.edu/~heister/ > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] error when trying to initialize PETSc preconditioners
Many thanks Timo. I tried PreconditionBlockJacobi and it worked. However, I am running a validation loop of progressively smaller meshes, and stops again with a error when I make the call to initialize for the second time (I have to do it, since the matrix is the same ... but has changed, I have to recompute the incomplete factorization), which is expected since ''initialize()'' can only be called once. A simple solution for this problem: declaring the preconditioner inside the validation loop rather than outside and pass it as an argument to the routine that solves the linear system. Yes, is there a better solution? Say for instance, as with trilinos preconditioners that have a ''clear()'' function member so that you can reinitialize it with a new matrix (this is particularly useful in adaptivity loops, like in step-32). That's all for now. Many Thanks. Ignacio. On Friday, August 12, 2016 at 9:31:36 PM UTC-5, Timo Heister wrote: > > Hey Ignacio, > > good to hear from you! > > Sadly this is not documented well, but the ICC preconditioner is only > implemented for BAIJ matrices and not MPIAIJ (which we use even if you > run on one processor), see > https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCICC.html > > You should be able to make it work by using a > PETScWrappers::SparseMatrix instead of an > PETScWrappers::MPI::SparseMatrix. Or, probably a better idea, just use > a preconditioner that works with parallel matrices like > PreconditionBlockJacobi. > > Best, > Timo > > > On Fri, Aug 12, 2016 at 9:46 PM, Ignacio Tomas > wrote: > > Hi, everybody. I am just trying to initialize either ICC, Jacobi or > > ParaSalis preconditioner from PETSc, using a PETSc matrix in parallel. > So > > far I am running just one mpi process. After the call to ''.initialize'' > I > > get: > > > > > > > > > > [0]PETSC ERROR: No support for this operation for this object type > > [0]PETSC ERROR: Matrix format mpiaij does not have a built-in PETSc ICC > > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html > for > > trouble shooting. > > [0]PETSC ERROR: Petsc Release Version 3.5.3, Jan, 31, 2015 > > [0]PETSC ERROR: ./main on a arch-linux2-c-opt named nachote by nachotex > Fri > > Aug 12 20:28:49 2016 > > [0]PETSC ERROR: Configure options --download-fblaslapack=1 > > --with-shared-libraries=1 --download-hypre=1 --download-mumps=1 > > --download-spooles=1 --download-scalapack=1 --download-metis=1 > > --download-parmetis=1 --download-blacs=1 --with-debugging=0 --with-x=0 > > [0]PETSC ERROR: #1 MatGetFactor() line 3960 in > > > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/mat/interface/matrix.c > > [0]PETSC ERROR: #2 PCSetup_ICC() line 18 in > > > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/ksp/pc/impls/factor/icc/icc.c > > > > [0]PETSC ERROR: #3 PCSetUp() line 902 in > > > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/ksp/pc/interface/precon.c > > > terminate called after throwing an instance of > > 'dealii::LACExceptions::ExcPETScError' > > what(): > > > > An error occurred in line <384> of file > > > > > > > in function > > void dealii::PETScWrappers::PreconditionICC::initialize(const > > dealii::PETScWrappers::MatrixBase&, const > > dealii::PETScWrappers::PreconditionICC::AdditionalData&) > > The violated condition was: > > ierr == 0 > > The name and call sequence of the exception was: > > ExcPETScError(ierr) > > Additional Information: > > An error with error number 56 occurred while calling a PETSc function > > > > > > > > > > Any clue? Did I forget to install something? > > > > -- > > The deal.II project is located at http://www.dealii.org/ > > For mailing list/forum options, see > > https://groups.google.com/d/forum/dealii?hl=en > > --- > > You received this message because you are subscribed to the Google > Groups > > "deal.II User Group" group. > > To unsubscribe from this group and stop receiving emails from it, send > an > > email to dealii+un...@googlegroups.com . > > For more options, visit https://groups.google.com/d/optout. > > > > -- > Timo Heister > http://www.math.clemson.edu/~heister/ > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Daniel Arndt and Jean-Paul Pelteret appointed deal.II Developers
All, I want to share the good news that Daniel Arndt and Jean-Paul Pelteret have accepted our invitation to become deal.II Developers. They are now listed at https://dealii.org/authors.html Both of them have provided patches throughout the library for several years already, but you all probably know them best from their tireless efforts in keeping the mailing list well organized and making sure that most of the questions get answers. Daniel & Jean-Paul, many thanks already for your past contributions, and we look forward to many more in the future! Happy hacking Wolfgang Bangerth -- Wolfgang Bangerth email:bange...@colostate.edu www: http://www.math.tamu.edu/~bangerth/ -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] error when trying to initialize PETSc preconditioners
Hey Ignacio, good to hear from you! Sadly this is not documented well, but the ICC preconditioner is only implemented for BAIJ matrices and not MPIAIJ (which we use even if you run on one processor), see https://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCICC.html You should be able to make it work by using a PETScWrappers::SparseMatrix instead of an PETScWrappers::MPI::SparseMatrix. Or, probably a better idea, just use a preconditioner that works with parallel matrices like PreconditionBlockJacobi. Best, Timo On Fri, Aug 12, 2016 at 9:46 PM, Ignacio Tomas wrote: > Hi, everybody. I am just trying to initialize either ICC, Jacobi or > ParaSalis preconditioner from PETSc, using a PETSc matrix in parallel. So > far I am running just one mpi process. After the call to ''.initialize'' I > get: > > > > > [0]PETSC ERROR: No support for this operation for this object type > [0]PETSC ERROR: Matrix format mpiaij does not have a built-in PETSc ICC > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for > trouble shooting. > [0]PETSC ERROR: Petsc Release Version 3.5.3, Jan, 31, 2015 > [0]PETSC ERROR: ./main on a arch-linux2-c-opt named nachote by nachotex Fri > Aug 12 20:28:49 2016 > [0]PETSC ERROR: Configure options --download-fblaslapack=1 > --with-shared-libraries=1 --download-hypre=1 --download-mumps=1 > --download-spooles=1 --download-scalapack=1 --download-metis=1 > --download-parmetis=1 --download-blacs=1 --with-debugging=0 --with-x=0 > [0]PETSC ERROR: #1 MatGetFactor() line 3960 in > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/mat/interface/matrix.c > [0]PETSC ERROR: #2 PCSetup_ICC() line 18 in > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/ksp/pc/impls/factor/icc/icc.c > [0]PETSC ERROR: #3 PCSetUp() line 902 in > /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/ksp/pc/interface/precon.c > terminate called after throwing an instance of > 'dealii::LACExceptions::ExcPETScError' > what(): > > An error occurred in line <384> of file > > in function > void dealii::PETScWrappers::PreconditionICC::initialize(const > dealii::PETScWrappers::MatrixBase&, const > dealii::PETScWrappers::PreconditionICC::AdditionalData&) > The violated condition was: > ierr == 0 > The name and call sequence of the exception was: > ExcPETScError(ierr) > Additional Information: > An error with error number 56 occurred while calling a PETSc function > > > > > Any clue? Did I forget to install something? > > -- > The deal.II project is located at http://www.dealii.org/ > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > --- > You received this message because you are subscribed to the Google Groups > "deal.II User Group" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to dealii+unsubscr...@googlegroups.com. > For more options, visit https://groups.google.com/d/optout. -- Timo Heister http://www.math.clemson.edu/~heister/ -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] error when trying to initialize PETSc preconditioners
Hi, everybody. I am just trying to initialize either ICC, Jacobi or ParaSalis preconditioner from PETSc, using a PETSc matrix in parallel. So far I am running just one mpi process. After the call to ''.initialize'' I get: [0]PETSC ERROR: No support for this operation for this object type [0]PETSC ERROR: Matrix format mpiaij does not have a built-in PETSc ICC [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting. [0]PETSC ERROR: Petsc Release Version 3.5.3, Jan, 31, 2015 [0]PETSC ERROR: ./main on a arch-linux2-c-opt named nachote by nachotex Fri Aug 12 20:28:49 2016 [0]PETSC ERROR: Configure options --download-fblaslapack=1 --with-shared-libraries=1 --download-hypre=1 --download-mumps=1 --download-spooles=1 --download-scalapack=1 --download-metis=1 --download-parmetis=1 --download-blacs=1 --with-debugging=0 --with-x=0 [0]PETSC ERROR: #1 MatGetFactor() line 3960 in /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/mat/interface/matrix.c [0]PETSC ERROR: #2 PCSetup_ICC() line 18 in /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/ksp/pc/impls/factor/icc/icc.c [0]PETSC ERROR: #3 PCSetUp() line 902 in /home/nachotex/Documents/deal/PETSC/petsc-3.5.3/src/ksp/pc/interface/precon.c terminate called after throwing an instance of 'dealii::LACExceptions::ExcPETScError' what(): An error occurred in line <384> of file in function void dealii::PETScWrappers::PreconditionICC::initialize(const dealii::PETScWrappers::MatrixBase&, const dealii::PETScWrappers::PreconditionICC::AdditionalData&) The violated condition was: ierr == 0 The name and call sequence of the exception was: ExcPETScError(ierr) Additional Information: An error with error number 56 occurred while calling a PETSc function Any clue? Did I forget to install something? -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Data exchanges in IO
Bruce, The p_local_solution vector is filled with values as the result of solving > a linear equation, p_local_rho is filled with values by looping over > individual elements and assigning values calculated from other quantities. > Part of the reason the p_local_rho vector doesn't have ghost nodes is that > you can't access individual elements for a ghosted vector. > You can't write to a ghosted vector, but you can read from it. Is this what you mean? The DataOut object appears to only handle ghosted vectors. It also appears > that you need to call compress on non-ghosted vectors (this was determined > empirically based on error messages). > Yes, for proper output you should use ghosted vectors that store the solution on all the dofs of locally owned cells. If you modify a parallel vector you have to call compress as mentioned in step-40 [1] for example. Would you want to have this information at some different place as well. If yes, where? > The vecScale function is just designed to multiply all values of an > un-ghosted vector by a constant. > You could also use operator*= [2] for this. If you don't change the values of the vectors using vecScale then the data > for "Phi" and "Rho" appears smooth (but the scale is unwieldy). If you > change the values of the vectors using vecScale, then "Phi" looks okay but > "Rho" appears to have unscaled values at locations that look like they are > at processor boundaries. I'm puzzled about how this is happening, since my > understanding is that a data exchange happens when you do the assignment > rho = p_local_rho. The vecScale function only operates on un-ghosted > vectors, but is it possible that it is missing some values somewhere? > This sounds weird as you do the exact same thing to (rho, p_local_rho) and (p_local_solution, phi). So you observe that p_local_solution looks wrong as well? Do you have a minimal working example showing this problem? Best, Daniel [1] https://dealii.org/8.4.0/doxygen/deal.II/step_40.html#LaplaceProblemassemble_system [2] https://www.dealii.org/8.4.0/doxygen/deal.II/classPETScWrappers_1_1VectorBase.html#a37900779c6049418c39bacc1d44f4260 -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Data exchanges in IO
Hi,
I'm trying to write out some data from a calculation being run in parallel.
I've got two vectors, p_local_solution and p_local_rho of type
PetscWrappers::MPI::Vector. The p_local_solution vector is initialized
using a statement
p_local_solution.reinit(p_locally_owned_dofs,p_locally_relevant_dofs,p_comm)
and should have ghost nodes included in it and the vector p_local_rho is
initialized as
p_local_rho.reinit(p_locally_owned_dofs,p_comm)
and doesn't have ghost nodes. The p_local_solution vector is filled with
values as the result of solving a linear equation, p_local_rho is filled
with values by looping over individual elements and assigning values
calculated from other quantities. Part of the reason the p_local_rho vector
doesn't have ghost nodes is that you can't access individual elements for
a ghosted vector. The code I'm using to write out the results looks like
the following:
DataOut data_out;
data_out.attach_dof_handler(p_dof);
LA::MPI::Vector phi;
phi.reinit(p_locally_owned_dofs,p_comm);
phi = p_local_solution;
vecScale(phi,1.0/p_units.EP_from_V);
phi.compress(VectorOperation::insert);
p_local_solution = phi;
data_out.add_data_vector(p_local_solution,"Phi");
LA::MPI::Vector rho;
vecScale(p_local_rho,1.0/p_charge);
rho.reinit(p_locally_owned_dofs,p_locally_relevant_dofs,p_comm);
p_local_rho.compress(VectorOperation::insert);
rho = p_local_rho;
data_out.add_data_vector(rho,"Rho");
The DataOut object appears to only handle ghosted vectors. It also appears
that you need to call compress on non-ghosted vectors (this was determined
empirically based on error messages). The vecScale function is just
designed to multiply all values of an un-ghosted vector by a constant. The
function has the form
void vecScale(LA::MPI::Vector &vec, double scale)
{
std::pair range;
range = vec.local_range();
unsigned int lo_idx = range.first;
unsigned int hi_idx = range.second;
for (unsigned int i=lo_idx; ihttp://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
Re: [deal.II] does the latest dealii prepackaged image file for Mac OSX come with 'other software packages' such as p4est, PETSc, Trilinos, etc?
Yes it does. Luca > On 11 ago 2016, at 20:36, thomas stephens wrote: > > It's not clear to me from the github dealii wiki for MacOSX whether or not > the prepackaged image file for OSX comes with all of the optional third party > libraries that are listed as optional software in the dealii README file. > > Tom > -- > The deal.II project is located at http://www.dealii.org/ > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > --- > You received this message because you are subscribed to the Google Groups > "deal.II User Group" group. > To unsubscribe from this group and stop receiving emails from it, send an > email to dealii+unsubscr...@googlegroups.com. > For more options, visit https://groups.google.com/d/optout. -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Vector-valued gradient of solution vector
Joel,
Yes the matrix you are assembling is a vector-valued mass matrix now.
For me your code is failing with
void dealii::FEValuesBase::get_function_gradients(const
InputVector&, std::vector >&) const [with InputVector =
dealii::Vector; int dim = 3; int spacedim = 3; typename
InputVector::value_type = double]
The violated condition was:
(fe->n_components()) == (1)
The name and call sequence of the exception was:
dealii::ExcDimensionMismatch((fe->n_components()),(1))
Additional Information:
Dimension 3 not equal to 1
and this is not surprising. What you are doing is to extract the gradients
of a vector-valued finite element solution. The object that should store
this should therefore be a std::vector> as you
want to store for each quadrature point and each component a Tensor<1,dim>.
What you want to do is really that in "Do not work". As you have two
DoFHandlers, you should also have two FEValues objects and two
cell_iterator corresponding to the correct DoFHandler. Then you can extract
the gradient of your scalar field and project it onto the ansatz space for
the vector-valued field. This should look like:
typename DoFHandler::active_cell_iterator cell_scalar =
dof_handler_scalar.begin_active();
typename DoFHandler::active_cell_iterator cell =
dof_handler.begin_active();
...
for (; cell!=endc; ++cell, ++cell_vector)
{
fe_values.reinit (cell);
fe_values_scalar.reinit (cell_scalar);
fe_values_scalar.get_function_gradients(test,fe_function_grad);
...
}
Best,
Daniel
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
[deal.II] Re: Vector-valued gradient of solution vector
Dear Daniel,
Thank you for a fast answer. I think I got the matrix and right-hand side
correct now, please see test.cc for the example code.
What I want to from the start was to get the vector-valued gradient from a
scalar field, but I could not get this to work. See section "Do not work"
in function setup_system in test.cc, this gives 0 everywhere. I want to use
the information from test vector but since it was made with a different
dof_handler it doesn´t seem to work. Is there a way around this problem?
Can a vector be transferred from one dof_handler to another? or is there a
better way to solve this?
Best,
Joel
On Friday, August 12, 2016 at 12:00:36 AM UTC+2, Daniel Arndt wrote:
>
> Joel,
>
> Does the equation you want to solve for have multiple components?
> Otherwise it would not make sense to multiply in the right hand side
> something with a gradient.
> If you want to project the gradient into a vector valued finite element
> ansatz space, then the matrix you are assembling looks weird.
> In what you are doing you are also considering the coupling between all
> components.
>
> You are calculating the gradient of your scalar finite element function
> correctly, but you need to find the correct component to use.
> For doing this you can use something like:
>
> const unsigned int component = fe.system_to_component_index(i).first;
> cell_rhs(i) += ((fe_values.shape_value (i, q_index) *
> fe_function_grad[q_index][component])*
> fe_values.JxW (q_index));
>
> If you want to assemble a mass matrix for a vector-valued finite element,
> you have also to restrict assembling for the matrix to the case
> where the component for the ith and jth ansatz function are the same.
>
> You might want to have a look at step-8 [1] and the module "Handling
> vector valued problems"[2].
>
> Best,
> Daniel
>
> [1]
> https://www.dealii.org/8.4.1/doxygen/deal.II/step_8.html#ElasticProblemassemble_system
> [2]
> https://www.dealii.org/8.4.1/doxygen/deal.II/group__vector__valued.html
>
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
/* -
*
* Copyright (C) 1999 - 2013 by the deal.II authors
*
* This file is part of the deal.II library.
*
* The deal.II library is free software; you can use it, redistribute
* it, and/or modify it under the terms of the GNU Lesser General
* Public License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
* The full text of the license can be found in the file LICENSE at
* the top level of the deal.II distribution.
*
* -
*
* Author: Wolfgang Bangerth, University of Heidelberg, 1999
*/
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
#include
using namespace dealii;
const double pi=3.14159265359;
double center_id=0;
template
class Step4
{
public:
Step4 ();
void run ();
private:
void make_grid ();
void setup_system();
void assemble_system ();
void solve ();
void output_results () const;
Triangulation triangulation;
FE_Qfe_scalar;
FESystem fe;
DoFHandler dof_handler_scalar;
DoFHandler dof_handler;
SparsityPattern sparsity_pattern;
SparseMatrix system_matrix;
ConstraintMatrix constraints;
Vector solution;
Vector system_rhs;
Vector test;
std::vector> nodes;
};
template
Step4::Step4 ()
:
fe_scalar (1),
fe (FE_Q(1),dim),
dof_handler (triangulation),
dof_handler_scalar (triangulation)
{}
template
void Step4::make_grid ()
{
GridGenerator::hyper_cube (triangulation, -2.5, 2.5);
triangulation.refine_global (4);
for (unsigned int step=0; step<2; ++step)
{
typename Triangulation::active_cell_iterator
cell = triangulation.begin_active(),
endc = triangulation.end();
for (; cell!=endc; ++cell)
{
for (unsigned int v=0;v < GeometryInfo::vertices_per_cell;++v)
{
Point center;
for (unsigned int i=0; ivertex(v));
if (std::fabs(distance_from_center) < 1.0)
{
cell->set_refine_flag ();
break;
}
}
}
triangulation.execute_coarsening_and_refinement ();
}
std::cout << " Number of active cells: "
<< triangulation.n_active_cells()
<< std::endl
<< " Total
[deal.II] Re: Extra layer around mesh
Ok Thanks, Joel On Thursday, August 11, 2016 at 6:05:41 PM UTC+2, Bruno Turcksin wrote: > > Joel, > > On Thursday, August 11, 2016 at 11:29:35 AM UTC-4, Joel Davidsson wrote: >> >> >> I just had a question about the third option. If one used FE_Nothing, >> what is the boundary that is used when doing calculations? Is it the >> boundary of the big box or the green box? Is this handle automatically or >> does one need to mark the green box with boundary indicator for example. >> > It will be the boundary on the big box and you can't use boundary > indicator because you are not on the boundary. The boundary is a geometric > "quantity", it doesn't care about the type of finite elements you are > using. So if you want to use Dirichlet "boundary" condition on the green > box, you need to build the ConstraintMatrix yourself. > > Best, > > Bruno > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. For more options, visit https://groups.google.com/d/optout.
