Square root of a function, or the associated vector, in C++

Question

I have computed a complex solution with the help of a mixed formulation space, i.e. I have a function $p$ in the space $W$, where

V = FiniteElement("Lagrange", triangle, 1)
W = V*V

The two components of $p$ represent real and imaginary part respectively, and I would now compute the modulus on my solution $$\left| p\right| = \sqrt{p_{real}^2 +p_{imag}^2}.$$

So far I have been able to compute $p_{real}^2 +p_{imag}^2$ via the code

Function modulus1(W->sub(0)->collapse());
Function modulus2(W->sub(0)->collapse());
Function modulus(W->sub(0)->collapse());

*modulus1.vector() =  *p_real.vector();
*modulus1.vector() *= *p_real.vector();
*modulus2.vector() =  *p_imag.vector();
*modulus2.vector() *= *p_imag.vector();

*modulus.vector() =  *modulus1.vector();
*modulus.vector() += *modulus2.vector();

but I still have to calculate the square root of modulus.

I tried accessing every element in modulus with a for loop, but it seems like the operator [] is not implemented.

Answer 1

Hi, for PETSc vector you could use VecPow, see here.

Comments

Hi, thank you for your answer!

That's what I am looking for, so I did a test with

*modulus.vector() = *p_real.vector();
PetscErrorCode VecPow(*modulus1.vector(), 2.0);

but I get a warning saying

warning: expression list treated as compound expression in initializer [-fpermissive]
PetscErrorCode VecPow(*modulus1.vector(), 2.0);

and there is no change in my vector modulus1.

---

So, how can I convert a std::shared_ptr to a dolfin::PETScVector?

Here I found a very old mail archive with the suggestion

//if x is a GenericVector, then

PETScVector* xx = dynamic_cast<PETScVector*>(&x);
if (xx == 0) {error("Not a PETSc Vector.");}

So I tried with

PETScVector* vv = dynamic_cast<PETScVector*>(&modulus);

but the warning I get is quite eloquent:

warning: dynamic_cast of ‘dolfin::Function modulus’ to ‘class dolfin::PETScVector*’ 
can never succeed

Even trying

PETScVector* vv = dynamic_cast<PETScVector*>(&modulus.vector());

leads to the error

error: cannot dynamic_cast ‘& dolfin::Function::vector()()’ (of type ‘class 
std::shared_ptr<dolfin::GenericVector>*’) to type ‘class dolfin::PETScVector*’ 
(source type is not polymorphic)

---

Try doing as shown here.

---

Yup, running

auto test = as_type<PETScVector>(*modulus.vector()).vec();
VecPow(test, 0.5);

succesfully changed the elements of modulus to their square root.

Thank you for your help.