Help with Vector<Number> template instantiations

[dragn...@gmail.com]

Hello,

I need to use Vector FullMatrix, BlockVector and BlockSparseMatrix with different datatype than what is already provided and instantiated (float, double...). I have the following problem when instantiating with a different type: the linking fails due to the multiple definitions of the following functions:
vector.cc:-1: error: multiple definition of `dealii::Vector<std::complex<float> >::operator=(std::complex<float>)'
vector.cc:-1: error: multiple definition of `dealii::Vector<int>::lp_norm(int) const'

These two functions are explicitly instantiated in vector.templates.h for float and int types. But, if I instantiate Vector<myNumber> in my code these two definitions will appear in both files (vector.cc and my source file) even though I use a different type.
Could you please let me know if there is a reason to instantiate them in the header file rather than in the source vector.cc file?
Thank you.

Dragan

[Wolfgang Bangerth]

Dragan,

> I need to use Vector FullMatrix, BlockVector and BlockSparseMatrix with
> different datatype than what is already provided and instantiated (float,
> double...). I have the following problem when instantiating with a different
> type: the linking fails due to the multiple definitions of the following
> functions:
> vector.cc:-1: error: multiple definition of
> `dealii::Vector<std::complex<float> >::operator=(std::complex<float>)'
> vector.cc:-1: error: multiple definition of `dealii::Vector<int>::lp_norm(int)
> const'
>
> These two functions are explicitly instantiated in vector.templates.h for
> float and int types. But, if I instantiate Vector<myNumber> in my code these
> two definitions will appear in both files (vector.cc and my source file) even
> though I use a different type.

Correct. Do you *explicitly* instantiate these classes in your code? If you
do, why? If they're already explicitly instantiated in the library, you do not
need to add explicit instantiations in your own code.

> Could you please let me know if there is a reason to instantiate them in the
> header file rather than in the source vector.cc file?

?? Can you point out where it is instantiated in the .h file? I think the only
file is in the .inst file, which is #included in the .cc file.

Best
Wolfgang

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[dragn...@gmail.com]

Hello Wolfgang,

Yes, I do explicitly instantiate these templates in my code. The datatype is different from what is already defined in deal.II. Something like:
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector.templates.h>
template class Vector<myNumberClass>;

But these two guys:

dealii::Vector<std::complex<float> >::operator=(std::complex<float>)

and

dealii::Vector<int>::lp_norm(int) const

are explicitly instantiated in vector.templates.h (lines 855 and 1476). So, if the file vector.templates.h is included in two source files there will be a linker error. And vector.templates.h is already included in lac/vector.cc, therefore the multiple definition error mentioning that file. The solution is to either put them into some source file in the library (i.e. in vector.cc) or mark them as inline. Probably better to place them into a source file.
What do you think?

Dragan

[Wolfgang Bangerth]

On 09/10/2016 10:33 AM, dragn...@gmail.com wrote:
> Hello Wolfgang,
>
> Yes, I do explicitly instantiate these templates in my code. The datatype is
> different from what is already defined in deal.II. Something like:
> #include <deal.II/lac/vector.h>
> #include <deal.II/lac/vector.templates.h>
> template class Vector<myNumberClass>;
>
> But these two guys:
> dealii::Vector<std::complex<float> >::operator=(std::complex<float>)
> and
> dealii::Vector<int>::lp_norm(int) const
> are explicitly instantiated in vector.templates.h (lines 855 and 1476).

Ah, I see. Good point. Yes, moving these two functions into the .cc file and
leaving a *declaration* of the explicit instantiation in the .h file is the
right solution.

Want to try writing a patch for this? We'd be happy to merge it!

[dragn...@gmail.com]

Yes, no problem. I see those files changed a lot since the 8.4.1 release. Anyway, I see the two spots are the same as in 8.4.1. I'll send you the patch.

Dragan

[Wolfgang Bangerth]

On 09/10/2016 11:22 AM, dragn...@gmail.com wrote:
> Yes, no problem. I see those files changed a lot since the 8.4.1 release.
> Anyway, I see the two spots are the same as in 8.4.1. I'll send you the patch.

Fantastic, thanks!
W.

[dragn...@gmail.com]

Hello Wolfgang,

Here is the patch. I added one more change in functions Vector<Number>::reinit:
   if (omit_zeroing_entries == false)
    *this = static_cast<Number>(0);
The line above fails if the Number type is not primitive. If it is a class it can't cast. So it could be:
   if (omit_zeroing_entries == false)
    *this = Number(0);
Therefore it will call a constructor or simply create a primitive type. I do not know if there are some ramifications of changing that line.
Please have a look.

Dragan

[Wolfgang Bangerth]

Dragan,

> Here is the patch. I added one more change in functions Vector<Number>::reinit:
> if (omit_zeroing_entries == false)
> *this = static_cast<Number>(0);
> The line above fails if the Number type is not primitive. If it is a class it
> can't cast. So it could be:
> if (omit_zeroing_entries == false)
> *this = Number(0);
> Therefore it will call a constructor or simply create a primitive type. I do
> not know if there are some ramifications of changing that line.
> Please have a look.

Yes, all of this looks correct. I think the *declaration* of these explicit
specializations should be moved to the bottom of the .h file, though.

Do you know how to use github? This would make sure that you get credit for
the patch, as you should! If you want to learn how to use it, take a look at
lecture 32.8 here:
http://www.math.colostate.edu/~bangerth/videos.html
If you want me to apply the patch, just say so, and that's ok with me as well.

Best
Wolfgang

[dragn...@gmail.com]

Hello Wolfgang,

It is ok, you can apply patch yourself.

I actually also had to specialize many template functions to work with the new type (those related to matrices and vectors, BlockMatrixBase::add, SparseMatrix::add and also local_apply_boundary_conditions etc.). There are some functions that do not work well with the user-defined types like Vector<>::allocate(), Vector<>::deallocate() and Vector<adouble>::operator =(). They cause seg. faults. I am not sure why, perhaps because of mem.align you use (Utilities::System::posix_memalign) and it is fine for primitive types. Anyway that is something internal for my case.

I am finishing what I started doing and will write you about an interesting application of deal.II. It is mostly useful for multi scale modelling in chemical engineering but is general in nature. It is coupling of deal.II with an equation-based simulator: basically, using deal.II only to assemble matrices/rhs (including non-linear FE cases), generating equations and then solving one or more of these systems together with other differential and algebraic equations in a large DAE system. All that is done from python (although implemented in c++). An example could be Lithium-ion batteries: electrodes are typically made out of a porous material composed of large numbers of solid particles, and particles are in a electrolyte. Particles (transport modelled using the FE method) are coupled at the larger electrode length scale via electrolyte transport (using the finite volume method). Two software are fully coupled and the boundary conditions can be set using the equations from the other software etc.

Dragan

[Bruno Turcksin]

Dragan,

On Wednesday, September 14, 2016 at 7:34:46 AM UTC-4, dragn...@gmail.com wrote:

I am finishing what I started doing and will write you about an interesting application of deal.II. It is mostly useful for multi scale modelling in chemical engineering but is general in nature. It is coupling of deal.II with an equation-based simulator: basically, using deal.II only to assemble matrices/rhs (including non-linear FE cases), generating equations and then solving one or more of these systems together with other differential and algebraic equations in a large DAE system. All that is done from python (although implemented in c++). An example could be Lithium-ion batteries: electrodes are typically made out of a porous material composed of large numbers of solid particles, and particles are in a electrolyte. Particles (transport modelled using the FE method) are coupled at the larger electrode length scale via electrolyte transport (using the finite volume method). Two software are fully coupled and the boundary conditions can be set using the equations from the other software etc.

Sounds great. Let us know when you have something working.

Best,

Bruno

[Wolfgang Bangerth]

Dragan,

> It is ok, you can apply patch yourself.

I've committed under your name (which hopefully I got right):
https://github.com/dealii/dealii/pull/3114

> I actually also had to specialize many template functions to work with the new
> type (those related to matrices and vectors, BlockMatrixBase::add,
> SparseMatrix::add and also local_apply_boundary_conditions etc.). There are
> some functions that do not work well with the user-defined types like
> Vector<>::allocate(), Vector<>::deallocate() and Vector<adouble>::operator
> =(). They cause seg. faults. I am not sure why, perhaps because of mem.align
> you use (Utilities::System::posix_memalign) and it is fine for primitive
> types. Anyway that is something internal for my case.

Hm, good question what is happening there. Feel free to propose individual
patches for each of these cases if you can identify an underlying reason.

> I am finishing what I started doing and will write you about an interesting
> application of deal.II. It is mostly useful for multi scale modelling in
> chemical engineering but is general in nature. It is coupling of deal.II with
> an equation-based simulator: basically, using deal.II only to assemble
> matrices/rhs (including non-linear FE cases), generating equations and then
> solving one or more of these systems together with other differential and
> algebraic equations in a large DAE system. All that is done from python
> (although implemented in c++). An example could be Lithium-ion batteries:
> electrodes are typically made out of a porous material composed of large
> numbers of solid particles, and particles are in a electrolyte. Particles
> (transport modelled using the FE method) are coupled at the larger electrode
> length scale via electrolyte transport (using the finite volume method). Two
> software are fully coupled and the boundary conditions can be set using the
> equations from the other software etc.

That would make for a fantastic code gallery project!

[dragn...@gmail.com]

Hello Wolfgang,

Ok, good. The commit is fine. I'll see if I can pinpoint the exact cause of segmentation faults.
One more question, in fact a small inconsistency that is of my interest: the function VectorTools::interpolate_boundary_values is not templated for the generic Number type but for double only (argument std::map<types::global_dof_index,double>& boundary_values). Is there a reason for that?

Dragan

[Wolfgang Bangerth]

On 09/15/2016 05:12 AM, dragn...@gmail.com wrote:
>
> One more question, in fact a small inconsistency that is of my interest:
> the function VectorTools::interpolate_boundary_values is not templated
> for the generic Number type but for double only (argument
> std::map<types::global_dof_index,double>& boundary_values). Is there a
> reason for that?

Which version of deal.II are you using? In the current development
version, this is properly templated so that the value type of the
Function object is the same as the value type of the boundary_values map.

[dragn...@gmail.com]

True, I've checked in the master and it is templated for the Function<dim,Number> too. I use 8.4.1.
It's fine then. Thanks.