Merge branch 'develop' of https://gitlab.inria.fr/mosaic/work-in-progress/bvpy into develop

src/bvpy/ivp.py

@@ -33,7 +33,7 @@ class Time():
         self.current += dt
 
 
-class Scheme():
+class FirstOrderScheme():
     def __init__(self, vform, dt=1):
         self._vform = None
         self.previous_solution = None
@@ -41,6 +41,17 @@ class Scheme():
 
         self.set_vform(vform)
 
+        self._type = None
+
+    @property
+    def type(self):
+        return self._type
+
+    @type.setter
+    def type(self, value):
+        assert value in ['implicit', 'explicit']
+        self._type = value
+
     def set_previous_solution(self, previous):
         self.previous_solution = previous
 
@@ -58,8 +69,11 @@ class Scheme():
     def ds(self):
         return self._vform.ds
 
+class ImplicitEuler(FirstOrderScheme):
+    def __init__(self, vform, dt=1):
+        super().__init__(vform, dt)
+        self.type = 'implicit'
 
-class ImplicitEuler(Scheme):
     def apply(self, u, v, sol):
         lhs, rhs = self._vform.get_form(u, v, sol)
         lhs = fe.Constant(self.dt)*lhs
@@ -68,15 +82,15 @@ class ImplicitEuler(Scheme):
         lhs += u*v*self.dx
         rhs += self.previous_solution*v*self.dx
 
+        # self._vform.lhs, self._vform.rhs = lhs, rhs
         return lhs, rhs
 
-
 class IVP(BVP):
-    def __init__(self, domain=None, scheme=None,  bc=None):
-        super().__init__(domain=domain, vform=scheme._vform, bc=bc)
+    def __init__(self, domain=None, vform=None, scheme=None, bc=None):
+        super().__init__(domain=domain, vform=vform, bc=bc)
 
         self.time = Time()
-        self.scheme = scheme
+        self.scheme = scheme(vform)
         self.previous_solution = self.solution.copy(deepcopy=True)
         self.scheme.set_previous_solution(self.previous_solution)
 
@@ -84,31 +98,41 @@ class IVP(BVP):
         self.solution.interpolate(solution)
         self.previous_solution.interpolate(solution)
 
-        lhs, rhs = self.scheme.apply(self.trialFunction, self.testFunction, self.solution)
-        self.vform.lhs = lhs
-        self.vform.rhs = rhs
+        # self.scheme.apply(self.trialFunction, self.testFunction, self.solution)
+        # self.vform.lhs = lhs
+        # self.vform.rhs = rhs
 
     def set_time(self, time):
         self.time.current = time
 
-    def solve(self, **kwargs):
+    def solve(self, lhs, rhs, **kwargs):
         if isinstance(self.vform.rhs, Zero):
-            self.solver = NonlinearSolver()
+            self.solver = NonLinearSolver()
+            jacobian = fe.derivative(lhs, self.solution, self.trialFunction)
+            self.solver.set_problem(self.solution, lhs,
+                                    jacobian, self.dirichletCondition,
+                                    self.pointSource)
         else:
             self.solver = LinearSolver()
-            self.solver.set_problem(self.solution, self.vform.lhs,
-                                    self.vform.rhs , self.dirichletCondition,
+            self.solver.set_problem(self.solution, lhs,
+                                    rhs , self.dirichletCondition,
                                     self.pointSource)
 
         self._set_solver_parameter(**kwargs)
 
         self.solver.solve()
 
-    def step(self):
-        self.time.step(self.scheme.dt)
-        self.solve()
-        self.previous_solution.assign(self.solution)
-
+    def step(self, dt):
+        self.scheme.dt = dt
+        lhs, rhs = self.scheme.apply(self.trialFunction, self.testFunction, self.solution)
+        if self.scheme.type == 'implicit':
+            self.time.step(self.scheme.dt)
+            self.solve(lhs, rhs)
+            self.previous_solution.assign(self.solution)
+        elif self.scheme.type == 'explicit':
+            self.solve(lhs, rhs)
+            self.time.step(self.scheme.dt)
+            self.previous_solution.assign(self.solution)
 
 
 if __name__ =='__main__':
@@ -119,16 +143,16 @@ if __name__ =='__main__':
     from bvpy.utils.visu import plot
     # bvp = BVP()
 
-    rectangle = Rectangle(x=-2, y=-2, dx=4, dy=4, resolution=0.05)
+    rectangle = Rectangle(x=-2, y=-2, length=4, width=4, resolution=0.05)
     # bvp.set_domain(domain)
 
     poisson = PoissonForm(source=0)
 
     diri = dirichlet(0, boundary='all')
 
-    implicit_euler = ImplicitEuler(poisson, dt=0.04)
+    # implicit_euler = ImplicitEuler(poisson, dt=0.04)
 
-    ivp = IVP(rectangle, implicit_euler, diri)
+    ivp = IVP(rectangle, poisson, ImplicitEuler, diri)
     sol_0 = create_expression('exp(-a*pow(x, 2) - a*pow(y, 2))', degree=2, a=5)
     ivp.set_initial_solution(sol_0)
     ivp.set_time(0)
@@ -136,16 +160,5 @@ if __name__ =='__main__':
     with fe.XDMFFile('solution.xdmf') as f:
         f.write(ivp.solution, ivp.time.current)
         for i in range(10):
-            ivp.step()
+            ivp.step(0.04)
             f.write(ivp.solution, ivp.time.current)
-
-
-    # plot(ivp.solution)
-    #
-    # plot(ivp.solution)
-    # ivp.step()
-    # plot(ivp.solution)
-    # ivp.step()
-    # plot(ivp.solution)
-    # ivp.step()
-    # plot(ivp.solution)

src/bvpy/templates/vforms/elasticity.py

@@ -163,9 +163,9 @@ class LinearElasticForm(ElasticForm):
 
     def construct_form(self, u, v, sol):
         self.lhs = fe.Constant(self._thickness)\
-                   * fe.inner(self._stress(u), self._strain(v)) * fe.dx
+                   * fe.inner(self._stress(u), self._strain(v)) * self.dx
 
-        self.rhs = fe.dot(self._source_term, v) * fe.dx(domain=v.function_space().mesh())\
+        self.rhs = fe.dot(self._source_term, v) * self.dx(domain=v.function_space().mesh())\
                    + self._neumann_term(v
 
 
@@ -262,8 +262,8 @@ class HyperElasticForm(ElasticForm):
             I1 = fe.tr(C) * J**(-2 / 3)
             W = (self._mu / 2) * (I1 - 3) + (K / 4) * ((J - 1)**2 + (fe.ln(J))**2)
 
-        self.lhs = self._thickness * W * fe.dx
-        self.lhs -= fe.dot(self._source_term, sol) * fe.dx
+        self.lhs = self._thickness * W * self.dx
+        self.lhs -= fe.dot(self._source_term, sol) * self.dx
         self.lhs -= self._neumann_term(v)
 
         self.lhs = fe.derivative(self.lhs, sol, v)

src/bvpy/templates/vforms/poisson.py

@@ -39,8 +39,8 @@ class PoissonForm(AbstractVform):
             self._source_term = source
 
     def construct_form(self, u, v, sol):
-        self.lhs = fe.inner(fe.grad(u), fe.grad(v))*fe.dx
-        self.rhs = self._source_term*v*fe.dx + self._neumann_term(v)
+        self.lhs = fe.inner(fe.grad(u), fe.grad(v))*self.dx
+        self.rhs = self._source_term*v*self.dx + self._neumann_term(v)
         # self.rhs = self._source_term*v*fe.dx + self._neuman*v*self.ds(0)