Python Function.function_space方法代码示例

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import function_space [as 别名]
 def plot(self, **kwargs):
     func = self._fefunc
     # fix a bug in the fenics plot function that appears when
     # the maximum difference between data values is very small
     # compared to the magnitude of the data
     values = func.vector().array()
     diff = max(values) - min(values)
     magnitude = max(abs(values))
     if diff < magnitude * 1e-8:
         logger.warning("PLOT: function values differ only by tiny amount -> plotting as constant")
         func = Function(func.function_space())
         func.vector()[:] = values[0]
     plot(func, **kwargs)

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import function_space [as 别名]
class OSI(Field):

    @classmethod
    def default_params(cls):
        params = Field.default_params()
        params.update(
            finalize=True,
            )
        return params

    def add_fields(self):
        params = self.params.copy_recursive()
        params["save"] = False
        params["plot"] = False
        #params["callback"] = False
        #params.pop("finalize")

        fields = []
        #fields.append(WSS(params=params))
        #return fields

        f = TimeIntegral("WSS", params=params, label="OSI")
        fields.append(f)
        fields.append(Magnitude(f, params=params))

        f = Magnitude("WSS", params=params)
        fields.append(f)
        fields.append(TimeIntegral(f, params=params, label="OSI"))

        #f = TimeIntegral("WSS", label="OSI")
        #fields.append(f)
        #fields.append(Magnitude(f))

        #f = Magnitude("WSS")
        #fields.append(f)
        #fields.append(TimeIntegral(f, label="OSI"))

        return fields

    def before_first_compute(self, get):
        tau = get("WSS")
        self.osi = Function(tau.sub(0).function_space().collapse())

    def compute(self, get):
        # Requires the fields Magnitude(TimeIntegral("WSS", label="OSI")) and
        # TimeIntegral(Magnitude("WSS"), label="OSI")
        #self.mag_ta_wss = get("Magnitude_TimeIntegral_WSS_OSI")
        #self.ta_mag_wss = get("TimeIntegral_Magnitude_WSS_OSI")
        self.mag_ta_wss = get("Magnitude_TimeIntegral_WSS-OSI")
        self.ta_mag_wss = get("TimeIntegral_Magnitude_WSS-OSI")

        if self.params.finalize:
            return None
        elif self.mag_ta_wss == None or self.ta_mag_wss == None:
            return None
        else:
            expr = conditional(self.ta_mag_wss < 1e-15,
                               0.0,
                               0.5 * (1.0 - self.mag_ta_wss / self.ta_mag_wss))
            self.osi.assign(project(expr, self.osi.function_space()))
            return self.osi

    def after_last_compute(self, get):
        self.mag_ta_wss = get("Magnitude_TimeIntegral_WSS-OSI")
        self.ta_mag_wss = get("TimeIntegral_Magnitude_WSS-OSI")
        #print self.name, " Calling after_last_compute"

        expr = conditional(self.ta_mag_wss < 1e-15,
                           0.0,
                           0.5 * (1.0 - self.mag_ta_wss / self.ta_mag_wss))
        self.osi.assign(project(expr, self.osi.function_space()))

        return self.osi

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import function_space [as 别名]
class ObjectiveFunctional(LinearOperator):
    """
    Provides data misfit, gradient and Hessian information for the data misfit
    part of a time-independent symmetric inverse problem.
    """
    __metaclass__ = abc.ABCMeta

    # Instantiation
    def __init__(self, V, Vm, bc, bcadj, \
    RHSinput=[], ObsOp=[], UD=[], Regul=[], Data=[], plot=False, \
    mycomm=None):
        # Define test, trial and all other functions
        self.trial = TrialFunction(V)
        self.test = TestFunction(V)
        self.mtrial = TrialFunction(Vm)
        self.mtest = TestFunction(Vm)
        self.rhs = Function(V)
        self.m = Function(Vm)
        self.mcopy = Function(Vm)
        self.srchdir = Function(Vm)
        self.delta_m = Function(Vm)
        self.MG = Function(Vm)
        self.Grad = Function(Vm)
        self.Gradnorm = 0.0
        self.lenm = len(self.m.vector().array())
        self.u = Function(V)
        self.ud = Function(V)
        self.diff = Function(V)
        self.p = Function(V)
        # Define weak forms to assemble A, C and E
        self._wkforma()
        self._wkformc()
        self._wkforme()
        # Store other info:
        self.ObsOp = ObsOp
        self.UD = UD
        self.reset()    # Initialize U, C and E to []
        self.Data = Data
        self.GN = 1.0   # GN = 0.0 => GN Hessian; = 1.0 => full Hessian
        # Operators and bc
        LinearOperator.__init__(self, self.delta_m.vector(), \
        self.delta_m.vector()) 
        self.bc = bc
        self.bcadj = bcadj
        self._assemble_solverM(Vm)
        self.assemble_A()
        self.assemble_RHS(RHSinput)
        self.Regul = Regul
        # Counters, tolerances and others
        self.nbPDEsolves = 0    # Updated when solve_A called
        self.nbfwdsolves = 0    # Counter for plots
        self.nbadjsolves = 0    # Counter for plots
        self._set_plots(plot)
        # MPI:
        self.mycomm = mycomm
        try:
            self.myrank = MPI.rank(self.mycomm)
        except:
            self.myrank = 0

    def copy(self):
        """Define a copy method"""
        V = self.trial.function_space()
        Vm = self.mtrial.function_space()
        newobj = self.__class__(V, Vm, self.bc, self.bcadj, [], self.ObsOp, \
        self.UD, self.Regul, self.Data, False)
        newobj.RHS = self.RHS
        newobj.update_m(self.m)
        return newobj

    def mult(self, mhat, y):
        """mult(self, mhat, y): do y = Hessian * mhat
        member self.GN sets full Hessian (=1.0) or GN Hessian (=0.0)"""
        N = self.Nbsrc # Number of sources
        y[:] = np.zeros(self.lenm)
        for C, E in zip(self.C, self.E):
            # Solve for uhat
            C.transpmult(mhat, self.rhs.vector())
            self.bcadj.apply(self.rhs.vector())
            self.solve_A(self.u.vector(), -self.rhs.vector())
            # Solve for phat
            E.transpmult(mhat, self.rhs.vector())
            Etmhat = self.rhs.vector().array()
            self.rhs.vector().axpy(1.0, self.ObsOp.incradj(self.u))
            self.bcadj.apply(self.rhs.vector())
            self.solve_A(self.p.vector(), -self.rhs.vector())
            # Compute Hessian*x:
            y.axpy(1.0/N, C * self.p.vector())
            y.axpy(self.GN/N, E * self.u.vector())
        y.axpy(1.0, self.Regul.hessian(mhat))

    # Getters
    def getm(self): return self.m
    def getmarray(self):    return self.m.vector().array()
    def getmcopyarray(self):    return self.mcopy.vector().array()
    def getVm(self):    return self.mtrial.function_space()
    def getMGarray(self):   return self.MG.vector().array()
    def getMGvec(self):   return self.MG.vector()
    def getGradarray(self):   return self.Grad.vector().array()
    def getGradnorm(self):  return self.Gradnorm
#.........这里部分代码省略.........

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import function_space [as 别名]
def testA():
    print "TEST A"
    quadrature_degree_old = parameters["form_compiler"]["quadrature_degree"]
      
    # setup problem
    mesh1 = UnitSquare(5, 5)
    V1 = FunctionSpace(mesh1, 'CG', 1)
    mesh2 = UnitSquare(10, 10)
    V2 = FunctionSpace(mesh2, 'CG', 1)
    u1 = TrialFunction(V1)
    v1 = TestFunction(V1)
    u2 = TrialFunction(V2)
    v2 = TestFunction(V2)
    f = Constant(1.0)

    # define boundary conditions
#    bc = DirichletBC(V, 0.0, "near(x[0], 0.0) || near(x[0], 1.0)")
    def u0_boundary(x, on_boundary):
        return (x[0] <= DOLFIN_EPS - 1.0 or (x[0] >= 0.0 - DOLFIN_EPS and x[1] < 1.0 - DOLFIN_EPS)) and on_boundary
    #    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
    u0 = Constant(0.0)
    bc1 = DirichletBC(V1, u0, u0_boundary)
    bc2 = DirichletBC(V2, u0, u0_boundary)

    # solution vector 
    u1_h = Function(V1)
    u2_h = Function(V2)

    # iterate expressions
    for run in range(2):
        if run == 1:
            u1_h.vector()[:] = np.ones(u1_h.function_space().dim())     # used for integration in form
        for j, ex in enumerate(EX):
            for qdegree in range(quadrature_degree):
                if qdegree == 0:
                    qdegree = -1
                parameters["form_compiler"]["quadrature_degree"] = qdegree
    
                # forms
    #            b1 = f * v1 * dx
    #            b2 = f * v2 * dx
    #            b1 = ex * v1 * dx
    #            b2 = ex * v2 * dx
                b1 = inner(nabla_grad(ex), nabla_grad(ex)) * v1 * dx
                b2 = inner(nabla_grad(ex), nabla_grad(ex)) * v2 * dx
                f1 = inner(nabla_grad(ex), nabla_grad(ex)) * u1_h * dx
                f2 = inner(nabla_grad(ex), nabla_grad(ex)) * u1_h * dx
    
                a1 = inner(nabla_grad(u1), nabla_grad(v1)) * dx
                a2 = inner(nabla_grad(u2), nabla_grad(v2)) * dx
    #            a1 = ex * inner(nabla_grad(u1), nabla_grad(v1)) * dx
    #            a2 = ex * inner(nabla_grad(u2), nabla_grad(v2)) * dx
    #            a1 = inner(nabla_grad(ex), nabla_grad(ex)) * inner(nabla_grad(u1), nabla_grad(v1)) * dx
    #            a2 = inner(nabla_grad(ex), nabla_grad(ex)) * inner(nabla_grad(u2), nabla_grad(v2)) * dx
    
                # compute solution
                if run == 0:
                    solve(a1 == b1, u1_h, bc1)
                    solve(a2 == b2, u2_h, bc2)
                    print "[V1] norms for quad=", qdegree, "ex", j, ":", norm(u1_h, 'L2'), norm(u1_h, 'H1')
                    print "[V2] norms for quad=", qdegree, "ex", j, ":", norm(u2_h, 'L2'), norm(u2_h, 'H1')
                                
                # compute norms
                if run == 1:
                    N1 = sqrt(assemble(f1))
                    N2 = sqrt(assemble(f2))
                    print "[V1] norm for quad=", qdegree, "ex", j, ":", N1
                    print "[V2] norm for quad=", qdegree, "ex", j, ":", N2
    
            print "-------------------------------------------------------------------------"
    parameters["form_compiler"]["quadrature_degree"] = quadrature_degree_old