How to project() a complex-valued function?

[Kev Se]

Hi!
I've recently swapped from another FEM package to deal.ii and love it so far. I recently managed so far to get a correct solution out for my complex-valued, nonlinear PDE. So that's good! What I am a bit confused about is the following: How do I use the VectorTools::project() type of function for a general, complex-valued function?

I've tried something like

/*****************/

typedef std::complex<double> dcomp;

template <int dim>

class Delta_of_x : public Function<dim, dcomp>

{

    public:

        virtual dcomp value(const Point<dim> &p,

        const unsigned int component = 0) const override;

};

template <int dim>

dcomp Delta_of_x<dim>::value(const Point<dim> &p,

const unsigned int component) const

{     dcomp returner;

       if (component == 0){

          returner= p[0]*p[0]; // dummy value, this can be more complicated

       }

     return returner;

}

/*********************************/

The above function works fine and as it should when used e.g. in the weak form.

After solving I want to output my solution, and my function Delta as well. So I tried to use:

/******************************/

AffineConstraints<dcomp> constraints;

constraints.close();

VectorTools::project(dof_handler, constraints, QGauss<dim>(degree + 2), Delta_of_x<dim>(), Delta_profile);

/*****************************************************/

This compiles just fine, but during linking I get the error message:

############################################

/path/to/myProject/dg_real.cc:690: error: undefined reference to 'void dealii::VectorTools::project<1, dealii::Vector<std::complex<double> > const, 1>(dealii::DoFHandler<1, 1> const&, dealii::AffineConstraints<dealii::Vector<std::complex<double> > const::value_type> const&, dealii::Quadrature<1> const&, dealii::Function<1, dealii::Vector<std::complex<double> > const::value_type> const&, dealii::Vector<std::complex<double> > const&, bool, dealii::Quadrature<(1)-(1)> const&, bool)'
/path/to/myProject/dg_real.cc:690: error: undefined reference to 'void dealii::VectorTools::project<2, dealii::Vector<std::complex<double> > const, 2>(dealii::DoFHandler<2, 2> const&, dealii::AffineConstraints<dealii::Vector<std::complex<double> > const::value_type> const&, dealii::Quadrature<2> const&, dealii::Function<2, dealii::Vector<std::complex<double> > const::value_type> const&, dealii::Vector<std::complex<double> > const&, bool, dealii::Quadrature<(2)-(1)> const&, bool)'
collect2: error: ld returned 1 exit status
make[2]: *** [CMakeFiles/dg_real.dir/build.make:127: dg_real] Error 1
make[1]: *** [CMakeFiles/Makefile2:90: CMakeFiles/dg_real.dir/all] Error 2
make: *** [Makefile:91: all] Error 2

############################################

What is going wrong here? I took the code for projecting from an example, and tried to adapt it to complex, but it seems like something is missing, or there is a mistake I cannot understand. Thank you for any input

[Wolfgang Bangerth]

On 10/26/22 12:57, Kev Se wrote:
>
> I've recently swapped from another FEM package to deal.ii and love it so
> far. I recently managed so far to get a correct solution out for my
> complex-valued, nonlinear PDE. So that's good! What I am a bit confused
> about is the following: How do I use the VectorTools::project() type of
> function for a general, complex-valued function?
>

Kev: It looks like the function you're trying to call is implemented but
not instantiated for complex arguments. This wouldn't be the first
function for which this is the case. Can you come up with as small a
testcase that illustrates the problem? The program doesn't have to do
anything useful other than showing that you can't link.

Separately, though: Have you made sure that you are working with the
latest version of deal.II? This may (or may not) be a bug that is
already fixed.

Best
W.

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

[Kev Se]

Dear Wolfgang, 

thank you for your swift reply! I am using Version 9.4.0, compiled from source just last week - I checked the github as well but it seems like there is no newer version available. 

I've tried to make a minimal working example and will attach it. The code now doesn't do anything expect for defining the class Delta_of_x that inherits from Function,  and its value() function, and then tries to project said function into an FE vector of nodal values that I want to save to disk. The current profile Delta(x) is in real, so I tried to change all datatypes to double but it didn't seem to help with the linking error (and I would want complex values for my actual problem, eventually). I'm starting to suspect that I'm doing something wrong with the project function, but I don't understand what if so. Either way, please find the .cc and a makefile attached. I  started from one of the examples at some point so there is probably some unnecessary includes and similar things left that I forgot to remove, I hope the problem is still understandable. 

Thank you for your help, best

Kev

[Kev Se]

Hi, 

after some tinkering - and considering the example/tutorial step-23, where the same project method was supposed to be used - I think I found a solution. The only difference to my code was that my output_results method, where the project() function was called, was declared as a const function. Removing that keyword removed the linking error. So in the class definition in the above mwe.cc, one changes the declaration from

void output_results(const unsigned int cycle) const;

to 

void output_results(const unsigned int cycle);

and identical replacements for the definition.  I guess one could have concluded that by thinking harder about what that const was doing there. Anyway, thanks for the help and keep up the great work!

Kev  

[Wolfgang Bangerth]

On 11/1/22 01:57, Kev Se wrote:
> *** Caution: EXTERNAL Sender ***

>
> Hi,
>
> after some tinkering - and considering the example/tutorial step-23,
> where the same project method was supposed to be used - I think I found
> a solution. The only difference to my code was that my output_results
> method, where the project() function was called, was declared as a

> /const/ function. Removing that keyword removed the linking error. So in

> the class definition in the above mwe.cc, one changes the declaration from
> void output_results(const unsigned int cycle) const;
> to
> void output_results(const unsigned int cycle);
>
> and identical replacements for the definition.  I guess one could have
> concluded that by thinking harder about what that const was doing there.
> Anyway, thanks for the help and keep up the great work!

Kev -- I still had this on my TODO list, but you already found the
answer. The error message the compiler gives is not useful because it
doesn't point out the error. The issue is that the function being called
has this signature:

template <int dim, typename VectorType, int spacedim>
void
project(const DoFHandler<dim, spacedim> & dof,
const AffineConstraints<typename VectorType::value_type>
&constraints,
const Quadrature<dim> &
quadrature,
const Function<spacedim, typename VectorType::value_type>
&function,
VectorType & vec,
const bool enforce_zero_boundary = false,
const Quadrature<dim - 1> &q_boundary = (...),
const bool project_to_boundary_first = false);

The output, `vec`, of course needs to be a writable vector, but since
you are in a `const` member function, in this line

VectorTools::project(
dof_handler, constraints, QGauss<dim>(degree + 2), DeltaX,
Delta_profile);

your `Delta_profile` vector is `const` as well. This makes no logical
sense, but for the compiler, it just means that it matches
VectorType = const dealii::Vector<std::complex<double> >
for which the *template* makes sense, but no implementation would make
sense, and so deal.II doesn't provide any instantiation either.


The problem with these sorts of things is that C++ does not (or better:
did not) provide convenient ways to say "this template can only be used
if VectorType is not a 'const' type". If we were using C++20, we could
have written this as

template <int dim, typename VectorType, int spacedim>
void project(...)
requires (std::is_const_v(VectorType) == false);

which tells the compiler to only consider the template whenever
VectorType is not a `const` type. In that case, you would have gotten a
compiler error, rather than a linker error, and the problem would have
been clearer. But we're not using C++20 yet, so this will have to wait
for a few years :-)

[Kev Se]

Dear Wolfgang, 

thank you for the elaborate and detailed response. I think I can understand the problem at hand. The "fix" of removing the "const" keyword indeed solved my problem, but I realized only later that it only works when have datatype "double" - for complex numbers I still got some error at linking, which should then not be related to non-access to a constant vector. I ended up running into similar problems with complex numbers when trying to use Petsc wrappers at a later stage, so I've know decided to switch to a coupled set of real equations instead. It is mentioned several times in the tutorials that this is the easier route, that seems to be the case for me.

Either way, your comments were very helpful in getting a clearer picture of the issue I ran into - time to extend my knowledge of (modern) C++ it would seem. 

Thanks again and best regards,

Kev