dolfin team mailing list archive

[Bartosz Sawicki <sawickib@xxxxxxxxxxxxx>]

I've written simple demo of curl-curl equation. It is based on 
eddy-currents phenomena in low conducting materials. So far, in the demo 
directory, there are no examples of this kind of PDE, where Nedelec 
elements are used.
If you found it useful, please consider placing it in the repository.

The patch which I'm sending, contains cpp and python code. Cpp works 
fine, but demo.py gives wrong results. I have no idea, whats wrong in 
the python code. Please, take a look, this is probably bug.

cheers,
BArtosz

diff -r 15ceb7ed6ff9 demo/pde/curl-curl/cpp/CurrentDensity.form
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demo/pde/curl-curl/cpp/CurrentDensity.form	Mon Mar 30 13:23:23 2009 -0600
@@ -0,0 +1,20 @@
+# Copyright (C) 2009 Bartosz Sawicki.
+# Licensed under the GNU LGPL Version 2.1.
+#
+# First added:  2009-03-30
+# Last changed: 2009-03-30
+#
+# The bilinear form a(v, u) and linear form L(v) for
+# curl operator.
+#
+# Compile this form with FFC: ffc -l dolfin CurrentDensity.form
+
+eL = VectorElement("CG", "tetrahedron", 1)
+eN = FiniteElement("Nedelec", "tetrahedron", 0)
+
+v = TestFunction(eL)
+u = TrialFunction(eL)
+T = Function(eN)
+
+a = dot(v,u)*dx
+L = dot(v,curl(T))*dx
diff -r 15ceb7ed6ff9 demo/pde/curl-curl/cpp/EddyCurrents.form
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demo/pde/curl-curl/cpp/EddyCurrents.form	Mon Mar 30 13:23:23 2009 -0600
@@ -0,0 +1,20 @@
+# Copyright (C) 2009 Bartosz Sawicki.
+# Licensed under the GNU LGPL Version 2.1.
+#
+# First added:  2009-03-30
+# Last changed: 2009-03-30
+#
+# The bilinear form a(v, u) and linear form L(v) for
+# curl-curl eddy currents equation.
+#
+# Compile this form with FFC: ffc -l dolfin EddyCurrents.form
+
+eN = FiniteElement("Nedelec", "tetrahedron", 0)
+eL = VectorElement("Lagrange", "tetrahedron", 1)
+
+v = TestFunction(eN)
+u = TrialFunction(eN)
+dBdt = Function(eL)
+
+a = dot(rot(v), rot(u))*dx 
+L = dot(v, dBdt)*dx 
diff -r 15ceb7ed6ff9 demo/pde/curl-curl/cpp/SConstruct
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demo/pde/curl-curl/cpp/SConstruct	Mon Mar 30 13:23:23 2009 -0600
@@ -0,0 +1,16 @@
+# This is an example of a simple SConstruct file for building
+# programs against DOLFIN. To build this demo, just type 'scons'.
+
+import os, commands
+
+# Get compiler from pkg-config
+compiler = commands.getoutput('pkg-config --variable=compiler dolfin')
+
+# Create a SCons Environment based on the main os environment
+env = Environment(ENV=os.environ, CXX=compiler)
+
+# Get compiler flags from pkg-config
+env.ParseConfig('pkg-config --cflags --libs dolfin')
+
+# Program name
+env.Program('demo', 'main.cpp')
diff -r 15ceb7ed6ff9 demo/pde/curl-curl/cpp/main.cpp
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demo/pde/curl-curl/cpp/main.cpp	Mon Mar 30 13:23:23 2009 -0600
@@ -0,0 +1,98 @@
+// Copyright (C) 2009 Bartosz Sawicki.
+// Licensed under the GNU LGPL Version 2.1.
+//
+// First added:  2009-03-30
+// Last changed: 2009-03-30
+//
+// Eddy currents phenomena in low conducting body can be 
+// described using electric vector potential and curl-curl operator:
+//    \nabla \times \nabla \times T = - \frac{\partial B}{\partial t}
+// Electric vector potential defined as:
+//    \nabla \times T = J
+//
+// Boundary condition 
+//    J_n = 0, 
+//    T_t=T_w=0, \frac{\partial T_n}{\partial n} = 0
+// which is naturaly fulfilled for zero Dirichlet BC with Nedelec (edge) 
+// elements.
+
+#include <dolfin.h>
+#include "EddyCurrents.h"
+#include "CurrentDensity.h"
+  
+using namespace dolfin;
+
+int main()
+{
+  
+  // Homogenous external magnetic field (dB/dt)
+  class Source : public Function
+  {
+  public:
+    void eval(real* values, const double* x) const
+    {
+      values[0] = 0.0;
+      values[1] = 0.0;
+      values[2] = 1.0;
+    }
+  };
+
+  // Zero Dirichlet BC
+  class Zero : public Function
+  {
+  public:
+    void eval(real* values, const double* x) const
+    {
+      values[0] = 0.0;
+      values[1] = 0.0;
+      values[2] = 0.0;
+    }
+  };
+
+  // Everywhere on external surface
+  class DirichletBoundary: public SubDomain
+  {
+    bool inside(const real* x, bool on_boundary) const
+    {
+      return on_boundary;
+    }
+  };
+
+  // Create demo mesh
+  UnitSphere mesh(10);
+
+  // Define functions
+  Source dBdt;
+  Zero zero;
+  Function T;
+  Function J;
+
+  // Define function space and boundary condition
+  EddyCurrentsFunctionSpace V(mesh);
+  DirichletBoundary boundary;
+  DirichletBC bc(V, zero, boundary);
+
+  // Use forms to define variational problem
+  EddyCurrentsBilinearForm a(V,V);
+  EddyCurrentsLinearForm L(V);
+  L.dBdt = dBdt;
+  VariationalProblem problem (a, L,  bc);
+  
+  // Solve problem using default solver
+  problem.solve(T);
+
+  // Define variational problem for current density (J)
+  CurrentDensityFunctionSpace V1(mesh);
+  CurrentDensityBilinearForm a1(V1,V1);
+  CurrentDensityLinearForm L1(V1);
+  L1.T = T;
+  VariationalProblem problem1(a1, L1);
+
+  // Solve problem using default solver
+  problem1.solve(J);
+
+  // Plot solution
+  plot(J);
+
+  return 0;
+}
diff -r 15ceb7ed6ff9 demo/pde/curl-curl/python/demo.py
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demo/pde/curl-curl/python/demo.py	Mon Mar 30 13:23:23 2009 -0600
@@ -0,0 +1,64 @@
+""" Eddy currents phenomena in low conducting body can be 
+described using electric vector potential and curl-curl operator:
+   \nabla \times \nabla \times T = - \frac{\partial B}{\partial t}
+Electric vector potential defined as:
+   \nabla \times T = J
+
+Boundary condition:
+   J_n = 0, 
+   T_t=T_w=0, \frac{\partial T_n}{\partial n} = 0
+which is naturaly fulfilled for zero Dirichlet BC with Nedelec (edge) 
+elements.
+"""
+
+__author__ = "Bartosz Sawicki (sawickib@xxxxxxxxxxxxx)"
+__date__ = "2009-03-30"
+__copyright__ = "Copyright (C) 2009 Bartosz Sawicki"
+__license__  = "GNU LGPL Version 2.1"
+
+from dolfin import *
+
+# Create mesh
+mesh = UnitSphere(10)
+
+# Define function spaces
+PN = FunctionSpace(mesh, "Nedelec", 0)
+P1 = VectorFunctionSpace(mesh, "CG", 1)
+
+# Define test and trial functions
+v0 = TestFunction(PN)
+u0 = TrialFunction(PN)
+v1 = TestFunction(P1)
+u1 = TrialFunction(P1)
+
+# Define functions
+dBdt = Function(P1, ("0.0", "0.0", "1.0"))
+zero = Function(P1, ("0.0", "0.0", "0.0"))
+T = Function(PN)
+J = Function(P1)
+
+# Dirichlet boundary
+class DirichletBoundary(SubDomain):
+    def inside(self, x, on_boundary):
+        return on_boundary
+
+# Boundary condition
+bc = DirichletBC(P1, zero, DirichletBoundary())
+
+# Eddy currents equation (using potential T) 
+Teqn = (dot(rot(v0), rot(u0))*dx, -dot(v0, dBdt)*dx)
+Tproblem = VariationalProblem( Teqn[0], Teqn[1], bc)
+T = Tproblem.solve()
+
+# Current density equation
+Jeqn = (dot(v1, u1)*dx, dot(v1, curl(T))*dx)
+Jproblem = VariationalProblem( Jeqn[0], Jeqn[1])
+J = Jproblem.solve()
+
+# Plot solution
+plot(J)
+
+# Hold plot
+interactive()
+
+