dolfin team mailing list archive

[Bartosz Sawicki <sawickib@xxxxxxxxxxxxx>]

Marie Rognes wrote:
> Bartosz Sawicki wrote:
> > > > 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.
> > >
> > > (1) There seems to be a minus sign discrepancy between the two demos 
> > > in the EddyCurrents part.
> >
> > You are right, there should be minus in both demos.
> >
> >
> > > (2) In the python demo, replace
> > >
> > >    bc = DirichletBC(P1, zero, DirichletBoundary())
> > >
> > > by    bc = DirichletBC(PN, zero, DirichletBoundary())
> >
> > Oh, good! How could I miss that. This is the bug which I was looking 
> > for. Now both examples works perfectly. Thank you very much.
> >
> > Do you want me to send improved patch? Or you can change this two 
> > characters?
>
> I have no pushing privileges for dolfin -- so I guess someone higher up 
> the hierarchy will get it done when they have time.
> But, it is probably easiest  if you fix the bugs and send the perfect 
> patch.

The perfect patch is in the attachment.

cheers,

diff -r 124a218df971 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	Thu Apr 02 11:34:05 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 124a218df971 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	Thu Apr 02 11:34:05 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 124a218df971 demo/pde/curl-curl/cpp/SConstruct
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/demo/pde/curl-curl/cpp/SConstruct	Thu Apr 02 11:34:05 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 124a218df971 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	Thu Apr 02 11:34:05 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 124a218df971 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	Thu Apr 02 11:34:05 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-04-02"
+__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(PN, 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()
+
+