translate python grad(u).trans code to C++

Kai · October 10, 2021, 10:48pm

Hi all :),
i try to translate the following piece of python code

self.u = self.fes.TrialFunction()
self.v = self.fes.TestFunction()
# StiffnessMatrix depending on material type:
if (self.material.type == "linear"):
   self.A = BilinearForm(self.fes)
   self.A += 2 * self.material.mu * InnerProduct(1 / 2 * (grad(self.u) + grad(self.u).trans),
                                                        1 / 2 * (grad(self.v) + grad(self.v).trans)) * dx

into corresponding C++ code so far i have

ngcomp::ProxyNode u = fes->GetTrialFunction();
ngcomp::ProxyNode v = fes->GetTestFunction();

Flags flags_bfa;
std::shared_ptr<T_BilinearFormSymmetric<double>> bfa = make_shared<T_BilinearFormSymmetric<double>> (fes, "a", flags_bfa);

bfa->AddIntegrator(make_shared<SymbolicBilinearFormIntegrator>(
                   material.mu * InnerProduct(
                    0.5f*(u->Deriv() + TransposeCF(u->Deriv())),
                    0.5f*(v->Deriv() + TransposeCF(v->Deriv()))),
                   VOL, VOL));

but the TransposeCF leads to

terminate called after throwing an instance of 'ngcore::Exception'
  what():  Transpose of non-matrix called
00:46:50: Das Programm ist abgestürzt.

So my question is how can i transpose v->Deriv() and u->Deriv() ?

Thank you for any help

joachim · October 11, 2021, 3:24am

Hi Kai,

can you show us how you created the finite element space fes ? Is it really vector-valued ?

For debugging, you can insert a

cout << *(u->Deriv()) << endl;

before the exception is thrown.

Joachim

Kai · October 11, 2021, 12:27pm

Hi Joachim,
thanks for quick response, i used the wrong FESpace. Now iam using VectorH1FESpace but the next problem follows : the SymbolicBilinearFormIntegrator does not work with Vector valued Coefficient functions …

hope this is ok to attach my new qustion here in this thread ?

joachim · October 11, 2021, 12:47pm

yes, at the end you need a scalar-valued function for the SymbolicBilinearFormIntegrator.

You get it, for example, by taking InnerProduct(grad(u), grad(v)).
grad(u) and grad(v) are matrix-valued.

Kai · October 11, 2021, 4:07pm

But i do already the InnerProduct like InnerProdukt(grad(u) + Transpose(grad(u)), grad(v) + Transpose(grad(v)))

bfa->AddIntegrator(make_shared<SymbolicBilinearFormIntegrator>(
                   material.mu * InnerProduct(
                    0.5f*(u->Deriv() + TransposeCF(u->Deriv())),
                    0.5f*(v->Deriv() + TransposeCF(v->Deriv()))),
                   VOL, VOL));

but this returns a coef innerproduct, fix size = 9, real … i do not see whats wrong here … mathematically this should work ?

joachim · October 11, 2021, 4:25pm

can you do a

auto func = material.mu * InnerProduct(
                    0.5f*(u->Deriv() + TransposeCF(u->Deriv())),
                    0.5f*(v->Deriv() + TransposeCF(v->Deriv())));
cout << *func << endl;

and send the output ?

If this does not help, you may send the complete example.

joachim · October 11, 2021, 4:25pm

can you do a

auto func = material.mu * InnerProduct(
                    0.5f*(u->Deriv() + TransposeCF(u->Deriv())),
                    0.5f*(v->Deriv() + TransposeCF(v->Deriv())));
cout << *func << endl;

and send the output ?

If this does not help, you may send the complete example.

Kai · October 11, 2021, 6:27pm

This lines give me the output :

2.59398e+07*(coef innerproduct, fix size = 9, real
  coef scale 0.5, real, dims = 3 x 3
    coef binary operation '+', real, dims = 3 x 3
      coef trial-function diffop = grad, real, dims = 3 x 3
      coef Matrix transpose, real, dims = 3 x 3
        coef trial-function diffop = grad, real, dims = 3 x 3
  coef scale 0.5, real, dims = 3 x 3
    coef binary operation '+', real, dims = 3 x 3
      coef test-function diffop = grad, real, dims = 3 x 3
      coef Matrix transpose, real, dims = 3 x 3
        coef test-function diffop = grad, real, dims = 3 x 3
)

Hope this say’s somthing to you …

EDIT :

the python line

print(InnerProduct(1 / 2 * (grad(self.u) + grad(self.u).trans), 1 / 2 * (grad(self.v) + grad(self.v).trans)))

gives me :

coef innerproduct, fix size = 9, real
  coef scale 0.5, real, dims = 3 x 3
    coef binary operation '+', real, dims = 3 x 3
      coef trial-function diffop = grad, real, dims = 3 x 3
      coef Matrix transpose, real, dims = 3 x 3
        coef trial-function diffop = grad, real, dims = 3 x 3
  coef scale 0.5, real, dims = 3 x 3
    coef binary operation '+', real, dims = 3 x 3
      coef test-function diffop = grad, real, dims = 3 x 3
      coef Matrix transpose, real, dims = 3 x 3
        coef test-function diffop = grad, real, dims = 3 x 3

so this seems to be very similar … but in python there is DifferentialSymbol dx doted with the result of the InnerProdukt … the dx seems to be defined as DifferentialSymbol(VOL) … should i dot with this in C++ ?

I can not find where it is defined …

joachim · October 19, 2021, 8:43am

I can’t see anything strange in the output, pls send me the complete example and I will have a look.

Kai · October 21, 2021, 11:05am

I have managed to implement the C++ code such that i now get same results on stdout with C++ and Python, but only for free boundary conditions. In Python i set dirichlet boundary conditions with

dirichlet = [4, 27, 9, 21, 28, 1, 33, 205, 277, 192, 3, 5, 18, 328]
fes = VectorH1(mesh, order=2, dirichlet=dirichlet, complex=True)

the integer list “dirichlet” i got with clicking on mesh faces in netgen ui. In C++ i do

Flags flags_fes;
flags_fes.SetFlag("order", 2);
flags_fes.SetFlag("complex", true);
flags_fes.SetFlag("dirichlet", "4, 27, 9, 21, 28, 1, 33, 205, 277, 192, 3, 5, 18, 328");
auto fes = make_shared<VectorH1FESpace>(ma, flags_fes);
BitArray dirichlet({4, 27, 9, 21, 28, 1, 33, 205, 277, 192, 3, 5, 18, 328});
fes->SetDirichletBoundaries(dirichlet);
fes->Update();
fes->FinalizeUpdate();
fes->UpdateDofTables();
fes->PrintReport(std::cout);

But the results are the same as if i do computation without setting boundary conditions.

So is there some thing special to call in C++ to apply boundary conditions ?

Kai · October 21, 2021, 1:34pm

Fixed !

    Flags flags_fes;
    flags_fes.SetFlag("order", 2);
    flags_fes.SetFlag("complex", true);
    //dirichlet(4, 27, 9, 21, 28, 1, 33, 205, 277, 192, 3, 5, 18, 328, 188, 154, 126);
    Array<double> dirbnd = {4, 27, 9, 21, 28, 1, 33, 205, 277, 192, 3, 5, 18, 328, 188, 154, 126};
    flags_fes.SetFlag("dirichlet", dirbnd);
    auto fes = make_shared<VectorH1FESpace>(ma, flags_fes);