Working non-linear heat conduction

demos/StationaryHeatTransfer.py

@@ -13,7 +13,7 @@ length = 30e-3
 width = 5.4e-3
 mesh = RectangleMesh(Point(0., 0.), Point(length, width), 100, 10)
 
-V = FunctionSpace(mesh, "CG", 2)
+V = FunctionSpace(mesh, "CG", 1)
 T = Function(V, name="Temperature")
 
 def left(x, on_boundary):
@@ -35,77 +35,17 @@ ds = Measure("ds", subdomain_data=facets)
 material = mf.MFrontNonlinearMaterial("../materials/src/libBehaviour.so",
                                       "StationaryHeatTransfer",
                                       hypothesis="plane_strain")
-problem = mf.MFrontNonlinearProblem(T, material)
+problem = mf.MFrontNonlinearProblem(T, material, quadrature_degree=0)
 problem.bc = bc
-print(problem.material.data_manager.K.shape)
+flux = mf.ThermalFlux(T)
+problem.define_form(flux)
 
-# definition of the stiffness matrix and residual using standard FEniCS operations
-a_Newton = (inner(grad(v), dot(as_tensor(dj_ddgT), grad(T_)))+
-            0*inner(grad(v), as_vector(dj_ddT))*T_)*dxm
-# a_Newton = inner(grad(v), dot(as_tensor(dj_ddgT), grad(T_)))*dxm
-res = -inner(grad(T_), j)*dxm
+T.interpolate(Constant(Tl))
 
-#problem.a = inner(self.strain_variation(self.du), dot(self.Ct, self.strain_variation(self.u_)))*measure*self.dx
-#problam.L = inner(self.strain_variation(self.u_), self.stress)*measure*self.dx l value of the deformation gradient
-#  using `MFront` conventions
-T0.assign(Constant(Tl))
-local_project(grad(T0), WgTs, gT)
-# copy the strain values to `MGIS`
-m.s0.gradients[:, :] = gT.vector().get_local().reshape((m.n, sj))
-m.s1.gradients[:, :] = gT.vector().get_local().reshape((m.n, sj))
-
-it = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator
-mgis_bv.integrate(m, it, 0, 0, m.n);
-bs = sj*sj+sj*1
-
-dj_ddgT.vector().set_local(np.lib.stride_tricks.as_strided(m.K,shape = (m.K.size // bs, sj, sj),
-                                                           strides = (m.K.strides[1] * bs,
-                                                                      m.K.strides[1] * sj,
-                                                                      m.K.strides[1])).flatten())
-dj_ddgT.vector().apply("insert")
-dj_ddT.vector().set_local(np.lib.stride_tricks.as_strided(m.K[sj*sj:],
-                                                          shape = (m.K.size // bs, sj),
-                                                          strides = (m.K.strides[1] * bs,
-                                                                     m.K.strides[1])).flatten())
-dj_ddT.vector().apply("insert")
-print(dj_ddT.vector().norm("l2"))
-
-
-T.assign(T0)
 problem.solve(T.vector())
 
-niter = 0
-tol = 1.e-8
-Nitermax = 200
-while nRes/nRes0 > tol and niter < Nitermax:
-    solve(A, dT.vector(), Res, "mumps")
-    # update the current estimate of the displacement at the end of the time step
-    T1.assign(T1+dT)
-    # compute the current estimate of the strain at the end of the
-    # time step using `MFront` conventions
-    local_project(grad(T1), WgTs, gT)
-    # copy the strain values to `MGIS`
-    m.s1.gradients[:, :] = gT.vector().get_local().reshape((m.n, sj))
-    # integrate the behaviour
-    it = mgis_bv.IntegrationType.IntegrationWithConsistentTangentOperator
-    # no time step here
-    mgis_bv.integrate(m, it, 0, 0, m.n);
-    # getting the stress and consistent tangent operator back to
-    # the FEniCS world.
-    j.vector().set_local(m.s1.thermodynamic_forces.flatten())
-    j.vector().apply("insert")
-    dj_ddgT.vector().set_local(np.lib.stride_tricks.as_strided(m.K,shape = (m.K.size // bs, sj, sj),
-                                                               strides = (m.K.strides[1] * bs,
-                                                                          m.K.strides[1] * sj,
-                                                                          m.K.strides[1])).flatten())
-    dj_ddgT.vector().apply("insert")
-    dj_ddT.vector().set_local(np.lib.stride_tricks.as_strided(m.K[sj*sj:],
-                                                              shape = (m.K.size // bs, sj),
-                                                              strides = (m.K.strides[1] * bs,
-                                                                         m.K.strides[1])).flatten())
-    dj_ddT.vector().apply("insert")
-    # solve Newton system
-    A, Res = assemble_system(a_Newton, res, bc0)
-    nRes = Res.norm("l2")
-    print("    Residual:", nRes)
-    niter += 1
+x = np.linspace(0, length)
+import matplotlib.pyplot as plt
+plt.plot(x, np.array([T(xi, width/2) for xi in x]))
+#plt.figure()
+#plot(project(flux.function[0], FunctionSpace(mesh, "DG",0)))
\ No newline at end of file

mfront_wrapper/gradient_flux.py

@@ -70,7 +70,7 @@ class Var(Gradient):
 
 
 class Flux:
-    def __init__(self, gradients, dual_gradient=None, shape=None, name=None):
+    def __init__(self, gradients, shape=None, dual_gradient=None, name=None):
         if type(gradients)==list:
             self.gradients = gradients
         else:

mfront_wrapper/nonlinear_problem.py

@@ -52,7 +52,7 @@ class MFrontNonlinearProblem(NonlinearProblem):
 
         # tangent bilinear form
         dg = self.flux.dual_gradient.variation(self.u_)
-        tangent_terms = [inner(dg, dot(t, g.variation(self.du)))
+        tangent_terms = [inner(dg, t*g.variation(self.du))
                          for (t, g) in zip(self.flux.tangents, self.flux.gradients)]
         self.a = sum(tangent_terms)*self.axi_coeff*self.dx
         # residual form (internal forces)
@@ -146,17 +146,13 @@ class MFrontNonlinearProblem(NonlinearProblem):
             K = self.material.data_manager.K
             buff = 0
             shapes = [ufl.shape(t) for t in self.flux.tangents]
-            ntot = sum([s[0]*s[1] if len(s)==2 else s[0] for s in shapes])
-            for t in self.flux.tangents:
-                s1, s2 = ufl.shape(t)
-                print(s1, s2)
-                t.vector().set_local(np.lib.stride_tricks.as_strided(K[buff:],
-                                                          shape = (K.size // ntot, s1, s2),
-                                                          strides = (K.strides[1] * ntot,
-                                                                     K.strides[1] * s1,
-                                                                     K.strides[1])).flatten())
-                t.vector().set_local(K.flatten())
-                buff += s1*s2
+            flattened_shapes = [s[0]*s[1] if len(s)==2 else s[0] for s in shapes]
+            for (i,t) in enumerate(self.flux.tangents):
+                s = flattened_shapes[i]
+                print(s, buff, K.shape, K[137,:])
+                t.vector().set_local(K[:,buff:buff+s].flatten())
+                #t.vector().set_local(K.flatten())
+                buff += s
             print(self.flux.tangents[0].vector().get_local())
             # sizes = self.material.get_state_variable_sizes()
             # for (s, i) in self.state_variables: