homog2d.py 21.5 KB
Newer Older
1 2 3 4
# coding: utf-8

import dolfin as fe
import logging
5
import numpy as np
6

7
from ho_homog import periodicity
Baptiste Durand's avatar
Baptiste Durand committed
8
from ho_homog.materials import epsilon, sigma, cross_energy
9 10 11 12 13
from ho_homog.toolbox_FEniCS import (
    DOLFIN_KRYLOV_METHODS,
    DOLFIN_LU_METHODS,
    local_project,
)
14

15
SOLVER_METHOD = "mumps"
16

17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33
U_NAMES = ("U1", "U2")
E_NAMES = ("E11", "E22", "E12")
EG_NAMES = ("EG111", "EG221", "EG121", "EG112", "EG222", "EG122")
EGG_NAMES = tuple(
    [n.replace("EG", "EGG") + "1" for n in EG_NAMES]
    + [n.replace("EG", "EGG") + "2" for n in EG_NAMES]
)

NAMES_MACRO_FIELDS = dict(
    U=U_NAMES,
    E=E_NAMES,
    EG=EG_NAMES,
    EGG=EGG_NAMES,
    EGbis=tuple([n.replace("EG", "EGbis") for n in EG_NAMES]),
    EGGbis=tuple([n.replace("EGG", "EGGbis") for n in EGG_NAMES]),
)

34

35 36
logging.getLogger("UFL").setLevel(logging.DEBUG)
logging.getLogger("FFC").setLevel(logging.DEBUG)
37
# https://fenicsproject.org/qa/3669/how-can-i-diable-log-and-debug-messages-from-ffc/
38
logger = logging.getLogger(__name__)
39 40 41


class Fenics2DHomogenization(object):
42
    """ #TODOC """
43

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
    def __init__(self, fenics_2d_rve, **kwargs):
        """[summary]

        Parameters
        ----------
        object : [type]
            [description]
        fenics_2d_rve : [type]
            [description]
        element : tuple or dict
            Type and degree of element for displacement FunctionSpace
            Ex: ('CG', 2) or {'family':'Lagrange', degree:2}
        solver : dict
            Choose the type of the solver, its method and the preconditioner.
            An up-to-date list of the available solvers and preconditioners
            can be obtained with dolfin.list_linear_solver_methods() and
            dolfin.list_krylov_solver_preconditioners().
61

62 63
        """
        self.rve = fenics_2d_rve
64
        self.topo_dim = topo_dim = fenics_2d_rve.dim
65 66 67 68 69 70 71 72
        try:
            bottom_left_corner = fenics_2d_rve.bottom_left_corner
        except AttributeError:
            logger.warning(
                "For the definition of the periodicity boundary conditions,"
                "the bottom left corner of the RVE is assumed to be on (0.,0.)"
            )
            bottom_left_corner = np.zeros(shape=(topo_dim,))
73
        self.pbc = periodicity.PeriodicDomain.pbc_dual_base(
74
            fenics_2d_rve.gen_vect, "XY", bottom_left_corner, topo_dim
75
        )
76

77
        solver = kwargs.pop("solver", {})
78
        # {'type': solver_type, 'method': solver_method, 'preconditioner': preconditioner}
79 80 81
        s_type = solver.pop("type", None)
        s_method = solver.pop("method", SOLVER_METHOD)
        s_precond = solver.pop("preconditioner", None)
82 83 84 85 86 87 88 89 90
        if s_type is None:
            if s_method in DOLFIN_KRYLOV_METHODS.keys():
                s_type = "Krylov"
            elif s_method in DOLFIN_LU_METHODS.keys():
                s_type = "LU"
            else:
                raise RuntimeError("The indicated solver method is unknown.")
        self._solver = dict(type=s_type, method=s_method)
        if s_precond:
91
            self._solver["preconditioner"] = s_precond
92

93
        element = kwargs.pop("element", ("Lagrange", 2))
94
        if isinstance(element, dict):
95
            element = (element["family"], element["degree"])
96 97 98 99 100 101 102 103
        self._element = element

        # * Function spaces
        cell = self.rve.mesh.ufl_cell()
        self.scalar_FE = fe.FiniteElement(element[0], cell, element[1])
        self.displ_FE = fe.VectorElement(element[0], cell, element[1])
        strain_deg = element[1] - 1 if element[1] >= 1 else 0
        strain_dim = int(topo_dim * (topo_dim + 1) / 2)
104
        self.strain_FE = fe.VectorElement("DG", cell, strain_deg, dim=strain_dim)
105
        # Espace fonctionel scalaire
106 107 108
        self.X = fe.FunctionSpace(
            self.rve.mesh, self.scalar_FE, constrained_domain=self.pbc
        )
109
        # Espace fonctionnel 3D : deformations, notations de Voigt
110
        self.W = fe.FunctionSpace(self.rve.mesh, self.strain_FE)
111
        # Espace fonctionel 2D pour les champs de deplacement
112
        # TODO : reprendre le Ve défini pour l'espace fonctionnel mixte. Par ex: V = FunctionSpace(mesh, Ve)
113 114 115
        self.V = fe.VectorFunctionSpace(
            self.rve.mesh, element[0], element[1], constrained_domain=self.pbc
        )
116

117 118
        # * Espace fonctionel mixte pour la résolution :
        # * 2D pour les champs + scalaire pour multiplicateur de Lagrange
119

120
        # "R" : Real element with one global degree of freedom
121
        self.real_FE = fe.VectorElement("R", cell, 0)
122 123
        self.M = fe.FunctionSpace(
            self.rve.mesh,
124
            fe.MixedElement([self.displ_FE, self.real_FE]),
125 126
            constrained_domain=self.pbc,
        )
127

128
        # Define variational problem
129 130 131
        self.v, self.lamb_ = fe.TestFunctions(self.M)
        self.u, self.lamb = fe.TrialFunctions(self.M)
        self.w = fe.Function(self.M)
132 133

        # bilinear form
134 135 136 137 138
        self.a = (
            fe.inner(sigma(self.rve.C_per, epsilon(self.u)), epsilon(self.v)) * fe.dx
            + fe.dot(self.lamb_, self.u) * fe.dx
            + fe.dot(self.lamb, self.v) * fe.dx
        )
139
        self.K = fe.assemble(self.a)
140 141 142 143 144
        if self._solver["type"] == "Krylov":
            self.solver = fe.KrylovSolver(self.K, self._solver["method"])
        elif self._solver["type"] == "LU":
            self.solver = fe.LUSolver(self.K, self._solver["method"])
            self.solver.parameters["symmetric"] = True
Baptiste Durand's avatar
Baptiste Durand committed
145
        try:
146
            self.solver.parameters.preconditioner = self._solver["preconditioner"]
Baptiste Durand's avatar
Baptiste Durand committed
147 148
        except KeyError:
            pass
149
        # fe.info(self.solver.parameters, True)
150

151 152 153 154 155 156 157
        self.localization = dict()
        # dictionary of localization field objects,
        # will be filled up when calling auxiliary problems (lazy evaluation)
        self.ConstitutiveTensors = dict()
        # dictionary of homogenized constitutive tensors,
        # will be filled up when calling auxiliary problems and
        # averaging localization fields.
158

159
    def homogenizationScheme(self, model):
160 161 162 163
        DictOfLocalizationsU = {}
        DictOfLocalizationsSigma = {}
        DictOfLocalizationsEpsilon = {}
        DictOfConstitutiveTensors = {}
164

165
        DictOfLocalizationsU["U"] = self.LocalizationU()["U"]
166

167 168
        if model == "E" or model == "EG" or model == "EGG":
            DictOfLocalizationsU["E"] = self.LocalizationE()["U"]
169

170 171
            DictOfLocalizationsSigma["E"] = self.LocalizationE()["Sigma"]
            DictOfLocalizationsEpsilon["E"] = self.LocalizationE()["Epsilon"]
172

173 174 175 176
            # TODO Renommer les champs enregistrés :
            #  for E_case, field in zip(E_NAMES, DictOfLocalizationsU["E"]):
            #   .rename(f"loc_{EG_case}_u", "")

177 178 179 180
        if model == "EG":
            DictOfLocalizationsU["EGbis"] = self.LocalizationEG()["U"]
            DictOfLocalizationsSigma["EGbis"] = self.LocalizationEGbis()["Sigma"]
            DictOfLocalizationsEpsilon["EGbis"] = self.LocalizationEGbis()["Epsilon"]
181

182 183 184 185
        if model == "EGG":
            DictOfLocalizationsU["EG"] = self.LocalizationEG()["U"]
            DictOfLocalizationsSigma["EG"] = self.LocalizationEG()["Sigma"]
            DictOfLocalizationsEpsilon["EG"] = self.LocalizationEG()["Epsilon"]
186

187 188 189
            DictOfLocalizationsU["EGGbis"] = self.LocalizationEGG()["U"]
            DictOfLocalizationsSigma["EGGbis"] = self.LocalizationEGGbis()["Sigma"]
            DictOfLocalizationsEpsilon["EGGbis"] = self.LocalizationEGGbis()["Epsilon"]
190

191 192 193 194 195 196
        # Prepare structure of DictOfConstitutiveTensors dictionary
        keys = list(DictOfLocalizationsSigma.keys())
        for k1 in keys:
            DictOfConstitutiveTensors[k1] = {}
            for k2 in keys:
                DictOfConstitutiveTensors[k1][k2] = None
197

198
        nk = len(keys)
199
        for i in range(nk):
200
            k1 = keys[i]
201
            for j in range(i, nk):
202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
                k2 = keys[j]
                sig_fields_k1_cases = DictOfLocalizationsSigma[k1]
                eps_fields_k2_cases = DictOfLocalizationsEpsilon[k2]
                C = self.CrossEnergy(k1, k2, sig_fields_k1_cases, eps_fields_k2_cases)
                DictOfConstitutiveTensors[k1][k2] = C
                DictOfConstitutiveTensors[k2][k1] = C.T

        # Rename localization fields :
        for k in DictOfLocalizationsU.keys():
            for n, loc_f in zip(NAMES_MACRO_FIELDS[k], DictOfLocalizationsU[k]):
                loc_f.rename(f"loc_{n}_u", "")

        for k in DictOfLocalizationsSigma.keys():
            names = NAMES_MACRO_FIELDS[k]
            for i in range(len(DictOfLocalizationsSigma[k])):
                n = names[i]
                loc_f = DictOfLocalizationsSigma[k][i]
                try:
                    loc_f.rename(f"loc_{n}_sig", "")
                except AttributeError:
                    DictOfLocalizationsSigma[k][i] = local_project(loc_f, self.W)
                    DictOfLocalizationsSigma[k][i].rename(f"loc_{n}_sig", "")
                loc_f = DictOfLocalizationsEpsilon[k][i]
                try:
                    loc_f.rename(f"loc_{n}_eps", "")
                except AttributeError:
                    DictOfLocalizationsEpsilon[k][i] = local_project(loc_f, self.W)
                    DictOfLocalizationsEpsilon[k][i].rename(f"loc_{n}_eps", "")
230

231 232 233 234 235 236
        return (
            DictOfLocalizationsU,
            DictOfLocalizationsSigma,
            DictOfLocalizationsEpsilon,
            DictOfConstitutiveTensors,
        )
237

238 239
    def LocalizationU(self):
        try:
240
            out = self.localization["U"]
241
        except KeyError:
242

243 244 245 246 247 248 249 250 251 252 253
            u = [
                fe.interpolate(fe.Constant((1.0, 0.0)), self.V),
                fe.interpolate(fe.Constant((0.0, 1.0)), self.V),
            ]
            s = [
                fe.interpolate(fe.Constant((0.0, 0.0, 0.0)), self.W),
                fe.interpolate(fe.Constant((0.0, 0.0, 0.0)), self.W),
            ]

            out = {"U": u, "Sigma": s, "Epsilon": s}
            self.localization["U"] = out
254
        return out
255

256
    def LocalizationE(self):
257
        """
258 259 260
        return E localization fields (as function of macro E)


261
        """
262
        try:
263
            out = self.localization["E"]
264
        except KeyError:
265
            f = [fe.Constant((0, 0)), fe.Constant((0, 0)), fe.Constant((0, 0))]
266
            Fload = [fe.interpolate(fo, self.V) for fo in f]
267 268 269 270 271
            epsilon0 = [
                fe.Constant((1, 0, 0)),
                fe.Constant((0, 1, 0)),
                fe.Constant((0, 0, 1)),
            ]
272
            Epsilon0 = [fe.interpolate(eps, self.W) for eps in epsilon0]
273
            u, sigm, eps = self.genericAuxiliaryProblem(Fload, Epsilon0)
274 275
            out = {"U": u, "Sigma": sigm, "Epsilon": eps}
            self.localization["E"] = out
276
        return out
277

278
    def LocalizationEbis(self):
279
        """
280
        return E stress/strain localization fields $u^U\diads\delta$ (as function of macro E)
281

282

283
        """
284
        try:
285
            out = self.localization["Ebis"]
286
        except KeyError:
287 288 289 290 291
            epsilon0 = [
                fe.Constant((1, 0, 0)),
                fe.Constant((0, 1, 0)),
                fe.Constant((0, 0, 1)),
            ]
292
            Epsilon0 = [fe.interpolate(eps, self.W) for eps in epsilon0]
293 294

            out = {
295 296 297 298
                "Sigma": [sigma(self.rve.C_per, Epsilon0[i]) for i in range(3)],
                "Epsilon": Epsilon0,
            }
            self.localization["Ebis"] = out
299 300 301
        return out

    def LocalizationEG(self):
302
        """
303 304
        return GradE localization fields (as function of macro GradE)

305
        """
306
        try:
307
            out = self.localization["EG"]
308
        except KeyError:
309
            self.LocalizationE()
310

311 312
            Fload = self.Sigma2Fload(self.localization["E"]["Sigma"])
            Epsilon0 = self.Displacement2Epsilon0(self.localization["E"]["U"])
313

314
            u, sigm, eps = self.genericAuxiliaryProblem(Fload, Epsilon0)
315 316
            out = {"U": u, "Sigma": sigm, "Epsilon": eps}
            self.localization["EG"] = out
317
        return out
318

319
    def LocalizationEGbis(self):
320
        """
321
        return GradE stress/strain localization fields for the EG model $u^E\diads\delta$ (as function of macro GradE)
322

323
        """
324
        try:
325
            out = self.localization["EGbis"]
326
        except KeyError:
327 328
            # h = self.lamination.total_thickness
            self.LocalizationE()
329

330
            Epsilon0 = self.Displacement2Epsilon0(self.localization["E"]["U"])
331 332

            out = {
333 334 335 336
                "Sigma": [sigma(self.rve.C_per, Epsilon0[i]) for i in range(6)],
                "Epsilon": Epsilon0,
            }
            self.localization["EGbis"] = out
337
        return out
338

339
    def LocalizationEGG(self):
340
        """
341
        return GradGradE localization fields (as function of macro GradGradE)
342

343 344 345 346
        >>> a = LaminatedComposite.paganosLaminate([1.0, 1.5, 1.0], [np.pi/3, -np.pi/3, np.pi/3])
        >>> Lam = Laminate2D(a)
        >>> plt=plotTransverseDistribution(Lam.LocalizationEGG()['U'])
        >>> show()
347
        """
348
        try:
349
            out = self.localization["EGG"]
350
        except KeyError:
351
            self.LocalizationEG()
352

353 354
            Fload = self.Sigma2Fload(self.localization["EG"]["Sigma"])
            Epsilon0 = self.Displacement2Epsilon0(self.localization["EG"]["U"])
355

356
            u, sigm, eps = self.genericAuxiliaryProblem(Fload, Epsilon0)
357

358 359
            out = {"U": u, "Sigma": sigm, "Epsilon": eps}
            self.localization["EGG"] = out
360
        return out
361

362
    def LocalizationEGGbis(self):
363
        """
364
        return GradGradE stress/strain localization fields for the EGG model $u^EG\diads\delta$ (as function of macro GradGradE)
365

366 367 368 369
        >>> a = LaminatedComposite.paganosLaminate([1.0, 1.5, 1.0], [np.pi/3, -np.pi/3, np.pi/3])
        >>> Lam = Laminate2D(a)
        >>> plt=plotTransverseDistribution(Lam.LocalizationEGGbis()['Sigma'])
        >>> show()
370
        """
371
        try:
372
            out = self.localization["EGGbis"]
373 374
        except KeyError:
            # h = self.lamination.total_thickness
375
            self.LocalizationEG()
376

377
            Epsilon0 = self.Displacement2Epsilon0(self.localization["EG"]["U"])
378 379

            out = {
380 381 382 383
                "Sigma": [sigma(self.rve.C_per, Epsilon0[i]) for i in range(12)],
                "Epsilon": Epsilon0,
            }
            self.localization["EGGbis"] = out
384
        return out
385

386
    def LocalizationSG(self):
387
        """
388 389
        return GradSigma localization fields (as function of macro GradSigma)

390
        """
391
        try:
392 393
            out = self.localization["SG"]
        except KeyError:
394
            self.LocalizationE()
395 396 397 398 399 400 401
            C = self.CrossEnergy(
                "E",
                "E",
                self.localization["E"]["Sigma"],
                self.localization["E"]["Epsilon"],
            )
            BEps = self.localization["E"]["Sigma"]
402 403
            n = len(BEps)
            BSigma = [[BEps[i][j] * C for j in range(n)] for i in range(n)]
404
            # Assembling the constitutive matrix
405 406 407 408
            BSigma_i = []
            for i in range(n):
                BSigma_j = []
                for j in range(n):
409
                    B = fe.Function(self.X)
410
                    for k in range(n):
411
                        print(C[k, j])
412
                        print(BEps[i].sub(k))
413 414 415
                        B = fe.interpolate(
                            B + BEps[i].sub(k) * fe.Constant(C[k, j]), self.X
                        )
416 417 418 419
                    BSigma_j = BSigma_j + [B]
                BSigma_i = BSigma_i + [BSigma_j]

            Fload = self.Sigma2Fload(BSigma)
420
            epsilon0 = [fe.Constant((0, 0, 0))] * 6
421
            Epsilon0 = [fe.interpolate(eps, self.W) for eps in epsilon0]
422

423
            u, sigm, eps = self.genericAuxiliaryProblem(Fload, Epsilon0)
424

425 426
            out = {"U": u, "Sigma": sigm, "Epsilon": eps}
            self.localization["SG"] = out
427
        return out
428

429
    def CrossEnergy(self, Key1, Key2, Loc1, Loc2):
Baptiste Durand's avatar
Baptiste Durand committed
430 431 432 433 434
        """
        Calcul l'energie croisee des tenseurs de localisation Loc1 et Loc2.
        Si Loc1 est un champ de contraintes alors Loc2 doit etre un champ de
        deformations (et inversement)
        """
435
        try:
436
            C = self.ConstitutiveTensors[Key1][Key2]
437
        except KeyError:
438
            logger.info(f"Compute cross energy {Key1} and {Key2}")
439
            if Key1 == Key2:
Baptiste Durand's avatar
Baptiste Durand committed
440 441 442 443
                len_1 = len(Loc1)
                C = np.zeros((len_1, len_1))
                for k in range(len_1):
                    for l in range(k, len_1):
444 445 446 447
                        C[k, l] = (
                            cross_energy(Loc1[k], Loc2[l], self.rve.mesh)
                            / self.rve.rve_area
                        )
448
                        C[l, k] = C[k, l]
449
            else:
Baptiste Durand's avatar
Baptiste Durand committed
450 451 452 453 454
                len_1 = len(Loc1)
                len_2 = len(Loc2)
                C = np.zeros((len_1, len_2))
                for k in range(len_1):
                    for l in range(len_2):
455 456 457 458
                        C[k, l] = (
                            cross_energy(Loc1[k], Loc2[l], self.rve.mesh)
                            / self.rve.rve_area
                        )
459 460
            try:
                self.ConstitutiveTensors[Key1][Key2] = C
461
            except KeyError:
462
                self.ConstitutiveTensors[Key1] = {}
463 464
                self.ConstitutiveTensors[Key1][Key2] = C

465 466
            try:
                self.ConstitutiveTensors[Key2][Key1] = C.T
467
            except KeyError:
468 469 470 471
                self.ConstitutiveTensors[Key2] = {}
                self.ConstitutiveTensors[Key2][Key1] = C.T
        return C

472
    def genericAuxiliaryProblem(self, Fload, Epsilon0):
473 474 475
        """
        This function compute the auxiliary problem of order N
        given the sources (of auxiliary problem N-1).
476

477 478
        It returns new localizations
        """
479 480 481 482 483

        U2 = []
        S2 = []
        E2 = []
        for i in range(len(Fload)):
Baptiste Durand's avatar
Baptiste Durand committed
484
            logger.info(f"Progression : load {i+1} / {len(Fload)}")
485 486 487 488
            L = (
                fe.dot(Fload[i], self.v)
                + fe.inner(-sigma(self.rve.C_per, Epsilon0[i]), epsilon(self.v))
            ) * fe.dx
489

490
            # u_s = fe.Function(self.V)
491 492 493
            res = fe.assemble(L)
            # self.bc.apply(res) #TODO à tester. Pas nécessaire pour le moment, la ligne # K,res = fe.assemble_system(self.a,L,self.bc) était commentée.
            # self.solver.solve(K, u_s.vector(), res) #* Previous method
494
            self.solver.solve(self.w.vector(), res)
495
            # * More info : https://fenicsproject.org/docs/dolfin/1.5.0/python/programmers-reference/cpp/la/PETScLUSolver.html
496

Baptiste Durand's avatar
Baptiste Durand committed
497 498 499 500 501 502 503 504 505 506 507 508 509
            u_only = fe.interpolate(self.w.sub(0), self.V)
            # ? Essayer de faire fe.assign(self.u_only, self.w.sub(0)) ?
            # ? Pour le moment u_only vit dans V et V n'est pas extrait
            # ? de l'espace fonctionnel mixte. Est-ce que cela marcherait si V
            # ? est extrait de M ?
            U2.append(u_only)
            eps = epsilon(u_only) + Epsilon0[i]
            E2.append(local_project(eps, self.W))
            S2.append(local_project(sigma(self.rve.C_per, eps), self.W))
            # TODO : Comparer les options :
            # TODO         e2 = fe.Function(self.W); e2.assign(self.RVE.epsilon(u_s) + Epsilon0[i])
            # TODO         interpolation
            # TODO         projection (locale)
510

511
        return U2, S2, E2
512 513 514 515

    def anyOrderAuxiliaryProblem(self, order=1):

        f = [fe.Constant((0, 0)), fe.Constant((0, 0)), fe.Constant((0, 0))]
516
        Fload = [fe.interpolate(fo, self.V) for fo in f]
517

518 519 520 521 522
        epsilon0 = [
            fe.Constant((1, 0, 0)),
            fe.Constant((0, 1, 0)),
            fe.Constant((0, 0, 1)),
        ]
523
        Epsilon0 = [fe.interpolate(eps, self.W) for eps in epsilon0]
524

525 526 527
        U = []
        Sigma = []
        Epsilon = []
528

529
        for i in range(order):
530
            logger.info(i)
531
            if i > 0:
532 533 534
                logger.info("compute load")
                Fload = self.Sigma2Fload(Sigma[i - 1])
                Epsilon0 = self.Displacement2Epsilon0(U[i - 1])
535

536
            u, sigm, eps = self.genericAuxiliaryProblem(Fload, Epsilon0)
537
            U = U + [u]
538 539
            Sigma = Sigma + [sigm]
            Epsilon = Epsilon + [eps]
540

541
        return U, Sigma, Epsilon
542 543

    def Sigma2Fload(self, S):
544 545
        """genere le chargement volumique associe au localisateur en
        contraintes du probleme precedent"""
546 547 548 549 550
        Fload = []
        S_avout = []
        # Oter la valeur moyenne
        for i in range(len(S)):
            Sj = []
551 552
            for j in range(int(self.topo_dim * (self.topo_dim + 1) / 2)):
                S_av = fe.assemble(S[i][j] * fe.dx) / self.rve.mat_area
553 554
                S_av = fe.interpolate(fe.Constant(S_av), self.X)
                Sj = Sj + [S[i][j] - S_av]
555
            S_avout = S_avout + [Sj]
556 557
        # Allocation du chargement
        for i in range(len(S)):
558
            Fload = Fload + [fe.as_vector((S_avout[i][0], S_avout[i][2]))]
559

560
        for i in range(len(S)):
561
            Fload = Fload + [fe.as_vector((S_avout[i][2], S_avout[i][1]))]
562

563
        return Fload
564

565
    def Displacement2Epsilon0(self, U):
566 567
        """Converti le localisateur en deplacement en champ de predeformation
        a appliquer au probleme auxiliaire suivant"""
568
        Epsilon0 = []
569

570
        for i in range(len(U)):
571 572 573 574 575 576 577 578 579
            Epsilon0 = Epsilon0 + [
                fe.as_vector(
                    (
                        U[i][0],
                        fe.interpolate(fe.Constant(0.0), self.X),
                        U[i][1] / fe.sqrt(2),
                    )
                )
            ]
580

581
        for i in range(len(U)):
582 583 584 585 586 587 588 589 590
            Epsilon0 = Epsilon0 + [
                fe.as_vector(
                    (
                        fe.interpolate(fe.Constant(0.0), self.X),
                        U[i][1],
                        U[i][0] / fe.sqrt(2),
                    )
                )
            ]
591

592
        return Epsilon0