dolfin team mailing list archive

[Patricia Brunner <question145791@xxxxxxxxxxxxxxxxxxxxx>]

New question #145791 on DOLFIN:
https://answers.launchpad.net/dolfin/+question/145791

I'm trying to port a working Python file to a C++ file. Unfortunately, the code, which works in Phyton, does not work in C++ and I don't know why. The result in Python is a vector with a norm around 1, whereas in C++ I get a vector with norm NaN.
I have attached both the Python and the C++ code below (which consists of two .ufl files and a .cpp file). Can anyone tell me what I have to change in the C++ code?

Thanks for your help,
Patricia

************************************************************************

Python code:

from dolfin import *

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

mesh = UnitSquare(64,64)

V0 = FunctionSpace(mesh, "DG", 0)
V1 = FunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 2)
W  = V1 * V2

z = Function(V0)
z.vector()[:] = 1.0

kappa = Function(V0)
mu = Function(V0)
kappa.vector()[:] = 1.0
mu.vector()[:] = 0.1

(sigma,u) = TrialFunction(W)
(v,tau)   = TestFunction(W)

a0 = (kappa*dot(grad(u),grad(tau)) + kappa*dot(grad(sigma),grad(v))- sigma*tau + mu*sigma*v + mu*u*tau)*dx +0.25*(u*tau+sigma*v)*ds
f0 = z*tau*dx

p = Function(W.sub(1))
y = Function(W.sub(0))
yp = Function(W)
p.vector()[:] = 0.0

A, F = assemble_system(a0,f0)
solve(A,yp.vector(),F)
(y,p) = yp.split()
print yp.vector().norm("l2")

plot(p,title="p")
interactive()

***************************************************

main.cpp:

#include <dolfin.h>
#include "Predual.h"
#include "PredualV0.h"
#include <vector>

using namespace dolfin;
using namespace std;

void conditional( Function& resVec, Function condVec, Constant alpha );

int main()
{

   // Create mesh
   UnitSquare mesh(64,64);

   // Create function spaces
   PredualV0::FunctionSpace V0(mesh);
   Predual::FunctionSpace W(mesh);

   // values of z (homogeneous lightning) in the three subdomains
   Function z(V0);
   z.vector() = 0.0;
   // values of kappa in the three subdomains
   Function kappa(V0);
   kappa.vector() = 1.0;
   // values of mu in the three subdomains
   Function mu(V0);
   mu.vector() = 0.1;

   Predual::BilinearForm a(W, W);
   Predual::LinearForm f(W);

   Function yp(W);
   Function& p = yp[1];
   
   Matrix A;
   Vector F;

   f.z = z;
   a.kappa = kappa;
   a.mu = mu;

   assemble_system(A, F, a, f);
   solve(A,yp.vector(),F);
   std::cout << yp.vector().norm("l2") << std::endl;

   plot(p,"p");
    
   return 0;
}

Predual.ufl:
V0 = FiniteElement("DG", "triangle", 0)
V1  = FiniteElement("CG", "triangle", 1)
V2 = FiniteElement("CG", "triangle", 2)
W = V1 * V2

(sigma, u) = TrialFunctions(W)
(tau, v)   = TestFunctions(W)

z = Coefficient(V0)
kappa = Coefficient(V0)
mu = Coefficient(V0)

a0 = (kappa*dot(grad(u),grad(tau)) + kappa*dot(grad(sigma),grad(v))- sigma*tau + mu*sigma*v + mu*u*tau)*dx + 0.25*(u*tau + sigma*v)*ds
f0 = z*tau*dx

a = a0
L = f0

PredualV0.ufl:
element = FiniteElement("DG", "triangle", 0)


-- 
You received this question notification because you are a member of
DOLFIN Team, which is an answer contact for DOLFIN.