[cig-commits] commit: updated MADDs-4 P1 forms to FFC 0.7.1

[Mercurial]

changeset:   111:93fed2bfd487
user:        Marc Spiegelman <mspieg at ldeo.columbia.edu>
date:        Fri Dec 11 13:33:52 2009 -0500
files:       MADDs-4/cpp/forms/SolitaryWave2D_n3m0_P1.ufl MADDs-4/cpp/forms/SolitaryWave3D_n3m0_P1.ufl
description:
updated MADDs-4 P1 forms to FFC 0.7.1


diff -r 3f22b9b7f49d -r 93fed2bfd487 MADDs-4/cpp/forms/SolitaryWave2D_n3m0_P1.ufl
--- a/MADDs-4/cpp/forms/SolitaryWave2D_n3m0_P1.ufl	Wed Dec 09 23:49:43 2009 -0500
+++ b/MADDs-4/cpp/forms/SolitaryWave2D_n3m0_P1.ufl	Fri Dec 11 13:33:52 2009 -0500
@@ -10,16 +10,19 @@
 #
 # with n=3, m=0
 #
-# Compile this form with FFC: ffc -l dolfin PressureEquation.ufl.
+# Compile this form with FFC: ffc -l dolfin -O SolitaryWave2D_n3m0_P1.ufl
 
 P1 = FiniteElement("Lagrange", triangle, 1)
-M = P1 + P1
+ME = P1 + P1
 
 
-(v,  u)  = TestFunctions(M)
-(dp, df) = TrialFunctions(M) # deltas for pressure and porosity
-(p, f)   = Functions(M)      # solution from last iteration
+v,q  = TestFunctions(ME)
+du   = TrialFunction(ME) # deltas for pressure and porosity
+u    = Function(ME)      # solution from last iteration
 
+# Split mixed functions
+dp, df = split(du) 
+p,  f  = split(u)
 
 dt = Constant(triangle)  # time step
 hdt = 0.5*dt # half time step
@@ -28,30 +31,24 @@ hsquared = Constant(triangle) # (h/delta
 hsquared = Constant(triangle) # (h/delta)^2
 
 # facet normal for forcing term
-n = VectorConstant(triangle)
+n = P2.cell().n 
 
 #Permeability 
-def K(f):
-	return f*f*f
+K = f*f*f
 
-def dKdf(f):
-	return 3*f*f
-
-#inverse bulk viscosity
-def Xi(f):
-	return 1
-
-def dXidf(f):
-	return 0
+#inverse bulk viscosity function (assume constant bulk viscosity)
+Xi = hsquared
 	
-
-
-a1 = (K(f)*inner(grad(v),grad(dp)) + hsquared*Xi(f)*v*dp)*dx + ((dKdf(f)*inner(grad(v),grad(p)) + hsquared*dXidf(f)*v*p - dKdf(f)*v.dx(1))*df)*dx + dKdf(f)*df*v*n[1]*ds 
-a2 = -u*hdt*hsquared*Xi(f)*dp*dx + u*(1 - hdt*dXidf(f)*hsquared*p)*df*dx 
-L1 = (K(f)*inner(grad(v),grad(p)) +hsquared*v*Xi(f)*p - K(f)*v.dx(1))*dx + v*K(f)*n[1]*ds 
-L2 = u*(f - hdt*hsquared*(Xi(f)*p))*dx 
+# calculate linear form for residual F(u)
+L1 = (K*inner(grad(v),grad(p)) + v*Xi*p - K*v.dx(1))*dx + v*K*n[1]*ds 
+L2 = q*(f - hdt*Xi*p)*dx 
  
-a = a1 + a2
 L = L1 + L2
 
+# bilinear form for weak form of the Jacobian J(u) (using functional differentiation of L)
 
+a = derivative(L, u, du)
+
+
+
+
diff -r 3f22b9b7f49d -r 93fed2bfd487 MADDs-4/cpp/forms/SolitaryWave3D_n3m0_P1.ufl
--- a/MADDs-4/cpp/forms/SolitaryWave3D_n3m0_P1.ufl	Wed Dec 09 23:49:43 2009 -0500
+++ b/MADDs-4/cpp/forms/SolitaryWave3D_n3m0_P1.ufl	Fri Dec 11 13:33:52 2009 -0500
@@ -10,15 +10,19 @@
 #
 # with n=3, m=0
 #
-# Compile this form with FFC: ffc -l dolfin SolitaryWave3D_n3m0_P1.ufl
+# Compile this form with FFC: ffc -l dolfin -O SolitaryWave3D_n3m0_P2.ufl
 
-P1 = FiniteElement("Lagrange", tetrahedron, 1)
-M = P1 + P1
+P1 = FiniteElement("Lagrange", "tetrahedron", 1)
+ME = P1 + P1
 
 
-(v,  u)  = TestFunctions(M)
-(dp, df) = TrialFunctions(M) # deltas for pressure and porosity
-(p, f)   = Functions(M)      # solution from last iteration
+v,q  = TestFunctions(ME)
+du   = TrialFunction(ME) # deltas for pressure and porosity
+u    = Function(ME)      # solution from last iteration
+
+# Split mixed functions
+dp, df = split(du) 
+p,  f  = split(u)
 
 
 dt = Constant(tetrahedron)  # time step
@@ -28,30 +32,24 @@ hsquared = Constant(tetrahedron) # (h/de
 hsquared = Constant(tetrahedron) # (h/delta)^2
 
 # facet normal for forcing term
-n = VectorConstant(tetrahedron)
+n = P1.cell().n # facet normal for forcing term
 
 #Permeability 
-def K(f):
-	return f*f*f
+K = f*f*f
 
-def dKdf(f):
-	return 3*f*f
-
-#inverse bulk viscosity
-def Xi(f):
-	return 1
-
-def dXidf(f):
-	return 0
+#inverse bulk viscosity function (assume constant bulk viscosity)
+Xi = hsquared
 	
-
-
-a1 = (K(f)*inner(grad(v),grad(dp)) + hsquared*Xi(f)*v*dp)*dx + ((dKdf(f)*inner(grad(v),grad(p)) + hsquared*dXidf(f)*v*p - dKdf(f)*v.dx(2))*df)*dx + dKdf(f)*df*v*n[2]*ds 
-a2 = -u*hdt*hsquared*Xi(f)*dp*dx + u*(1. - hdt*dXidf(f)*hsquared*p)*df*dx 
-L1 = (K(f)*inner(grad(v),grad(p)) +hsquared*v*Xi(f)*p - K(f)*v.dx(2))*dx + v*K(f)*n[2]*ds 
-L2 = u*(f - hdt*hsquared*(Xi(f)*p))*dx 
+# calculate linear form for residual F(u)
+L1 = (K*inner(grad(v),grad(p)) + v*Xi*p - K*v.dx(2))*dx + v*K*n[2]*ds 
+L2 = q*(f - hdt*Xi*p)*dx 
  
-a = a1 + a2
 L = L1 + L2
 
+# bilinear form for weak form of the Jacobian J(u) (using functional differentiation of L)
 
+a = derivative(L, u, du)
+
+
+
+