Python LinearPDE.setValue方法代码示例

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
class SimpleStokesProblem(SaddlePointProblem):
      """
      simple example of saddle point problem
      """
      def __init__(self,domain):
         super(SimpleStokesProblem, self).__init__(self)

         self.__pde_u=LinearPDE(domain)
         self.__pde_u.setSymmetryOn()
         self.__pde_u.setValue(A=identityTensor4(dom))

         self.__pde_p=LinearPDE(domain)
         self.__pde_p.setReducedOrderOn()
         self.__pde_p.setSymmetryOn()
         self.__pde_p.setValue(D=1.)

      def initialize(self,f=Data(),fixed_u_mask=Data()):
         self.__pde_u.setValue(q=fixed_u_mask,Y=f)
      def inner(self,p0,p1):
         return integrate(p0*p1,Function(self.__pde_p.getDomain()))

      def solve_f(self,u,p,tol=1.e-8):
         self.__pde_u.setTolerance(tol)
         self.__pde_u.setValue(X=grad(u)+p*kronecker(self.__pde_u.getDomain()))
         return  self.__pde_u.getSolution()
      def solve_g(self,u,tol=1.e-8):
         self.__pde_p.setTolerance(tol)
         self.__pde_p.setValue(X=-u) 
         dp=self.__pde_p.getSolution()
         return  dp

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
    def setUpPDE(self):
        """
        Creates and returns the underlying PDE.

        :rtype: `LinearPDE`
        """
        if self.__pde is None:
            if not HAVE_DIRECT:
                raise ValueError("Either this build of escript or the current MPI configuration does not support direct solvers.")
            pde=LinearPDE(self.__domain, numEquations=2)
            D=pde.createCoefficient('D')
            A=pde.createCoefficient('A')
            A[0,:,0,:]=kronecker(self.__domain.getDim())
            A[1,:,1,:]=kronecker(self.__domain.getDim())
            pde.setValue(A=A, D=D)
            if self.__fixAtBottom:
                DIM=self.__domain.getDim()
                z = self.__domain.getX()[DIM-1]
                pde.setValue(q=whereZero(z-self.__BX[DIM-1][0])*[1,1])

            pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
            pde.getSolverOptions().setTolerance(self.__tol)
            pde.setSymmetryOff()
        else:
            pde=self.__pde
            pde.resetRightHandSideCoefficients()
        return pde

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
   def test_Poisson2Classic(self):
        x=Solution(self.domain).getX()
        # --- set exact solution ----
        u_ex=Data(1.,(2,),Solution(self.domain))
        g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
        A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
        for i in range(self.domain.getDim()):
           u_ex[0]+= 1*(i+1) *x[i]
           g_ex[0,i]=1*(i+1)
           u_ex[1]+= 2*(i+1)*x[i]
           g_ex[1,i]=2*(i+1)
           A[0,i,0,i]=1
           A[1,i,1,i]=1

        # create PDE:
        pde=LinearPDE(self.domain,numEquations=2)
        pde.setSymmetryOn()
        mask=whereZero(x[0])*[1,1]
        pde.setValue(r=u_ex,q=mask)
        pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
        # -------- get the solution ---------------------------
        pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
        pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
        pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
        pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
        if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
        if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
        if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
        pde.getSolverOptions().setAMGInterpolation(SolverOptions.CLASSIC_INTERPOLATION)

        u=pde.getSolution()
        # -------- test the solution ---------------------------
        error=Lsup(u-u_ex)/Lsup(u_ex)
        self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
   def test_PoissonWithDirectInterpolation(self):
        x=Solution(self.domain).getX()
        # --- set exact solution ----
        u_ex=Scalar(1,Solution(self.domain))
        g_ex=Vector(0.,Solution(self.domain))
        for i in range(self.domain.getDim()):
           u_ex+=(i+1)*x[i]
           g_ex[i]=(i+1)

        # create PDE:
        pde=LinearPDE(self.domain,numEquations=1)
        pde.setSymmetryOn()
        mask=whereZero(x[0])
        pde.setValue(r=u_ex,q=mask)
        pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
        # -------- get the solution ---------------------------
        pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
        pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
        pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
        pde.getSolverOptions().setNumPreSweeps(3)
        pde.getSolverOptions().setNumPostSweeps(3)
        pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
        if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
        if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
        if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
        pde.getSolverOptions().setAMGInterpolation(SolverOptions.DIRECT_INTERPOLATION)

        u=pde.getSolution()
        # -------- test the solution ---------------------------
        error=Lsup(u-u_ex)/Lsup(u_ex)
        self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
    def setUpPDE(self):
        """
        Return the underlying PDE.

        :rtype: `LinearPDE`
        """
        if self.__pde is None:
            dom=self.__domain
            x = dom.getX()
            DIM=dom.getDim()
            q=whereZero(x[DIM-1]-inf(x[DIM-1]))
            for i in range(DIM-1):
                xi=x[i]
                q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))

            pde=LinearPDE(dom, numEquations=1)
            pde.getSolverOptions().setTolerance(self.__tol)
            pde.setSymmetryOn()
            A=pde.createCoefficient('A')
            X=pde.createCoefficient('X')
            pde.setValue(A=A, X=X, q=q)

        else:
            pde=self.__pde
            pde.resetRightHandSideCoefficients()
        return pde

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
    def getPotential(self):
        """
        returns a list containing 3 lists one for each the primary, secondary
        and total potential.
        """


        primCon=self.primaryConductivity
        coords=self.domain.getX()
        pde=LinearPDE(self.domain, numEquations=1)
        tol=1e-8
        pde.getSolverOptions().setTolerance(tol)
        pde.setSymmetryOn()

        DIM=self.domain.getDim()
        x=self.domain.getX()
        q=whereZero(x[DIM-1]-inf(x[DIM-1]))
        for i in xrange(DIM-1):
            xi=x[i]
            q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
        A = self.secondaryConductivity * kronecker(self.domain)
        pde.setValue(A=A,q=q)

        delPhiSecondary = []
        delPhiPrimary = []
        delPhiTotal = []
        if(len(self.electrodes[0])==3):

            for i in range(self.numElectrodes-1):
                analyticRs=Data(0,(3,),ContinuousFunction(self.domain))
                analyticRs[0]=(coords[0]-self.electrodes[i][0])
                analyticRs[1]=(coords[1]-self.electrodes[i][1])
                analyticRs[2]=(coords[2])
                rsMag=(analyticRs[0]**2+analyticRs[1]**2+analyticRs[2]**2)**0.5
                analyticPrimaryPot=(self.current*(1./primCon))/(2*pi*(rsMag+(whereZero(rsMag)*0.0000001))) #the magic number 0.0000001 is to avoid devide by 0
                analyticRsPolePower=(analyticRs[0]**2+analyticRs[1]**2+analyticRs[2]**2)**1.5
                analyticRsPolePower = analyticRsPolePower+(whereZero(analyticRsPolePower)*0.0000001)
                gradUPrimary = Data(0,(3,),ContinuousFunction(self.domain))
                gradUPrimary[0] =(self.current/(2*pi*primCon)) * (analyticRs[0]/analyticRsPolePower)
                gradUPrimary[1] =(self.current/(2*pi*primCon)) * (analyticRs[1]/analyticRsPolePower)
                gradUPrimary[2] =(self.current/(2*pi*primCon)) * (analyticRs[2]/analyticRsPolePower)
                gradUPrimary=-gradUPrimary
                X=(primCon-self.secondaryConductivity) * gradUPrimary
                pde.setValue(X=X)
                u=pde.getSolution()
                loc=Locator(self.domain,self.electrodes[i+1])
                delPhiSecondary.append(loc.getValue(u))
                delPhiPrimary.append(loc.getValue(analyticPrimaryPot))
        else:
            raise NotImplementedError("2d forward model is not yet implemented")

        self.delPhiSecondary = delPhiSecondary
        self.delPhiPrimary = delPhiPrimary
        for i in range(len(delPhiPrimary)):
            delPhiTotal.append(delPhiPrimary[i] + delPhiSecondary[i])
        self.delPhiTotal=delPhiTotal
        return [delPhiPrimary, delPhiSecondary, delPhiTotal]

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
 def getMaxEigenvalue(self):
    """
    return the max eigenvalue of the model
    type: float
    """
    dom = self.getDomain()
    dim = dom.getDim()
    pdeKu_P = LinearPDE(dom,numEquations=dim,numSolutions=dim)
    T = self.getCurrentTangent()
    pdeKu_P.setValue(A=T)
    pdeKu_P.getOperator().saveMM('tempMTX.mtx')
    mtx = mmread('tempMTX.mtx')
    maxEigVal = max(eigs(mtx,k=1,return_eigenvectors=False,which='LR'))
    return maxEigVal.real

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
    def test_PDE2D(self):
        dx_tests=0.1
        sigma0=1.
        electrodes=[(0.5-2*dx_tests,1.), (0.5-dx_tests,1.), (0.5+dx_tests,1.), (0.5+2*dx_tests,1.)]
        domain=finRectangle(20,20, d1=mpisize,  diracPoints=electrodes, diracTags=["sl0", "sl1", "sr0", "sr1"] )
        loc=Locator(domain,electrodes[2:])

        # this creates some reference Data:
        x=domain.getX()
        q=whereZero(x[0]-inf(x[0]))+whereZero(x[0]-sup(x[0]))+whereZero(x[1]-inf(x[1]))
        ppde=LinearPDE(domain, numEquations=1)
        s=Scalar(0.,DiracDeltaFunctions(domain))
        s.setTaggedValue("sl0" ,1.)
        s.setTaggedValue("sl1",-1.)
        ppde.setValue(A=kronecker(2)*sigma0, q=q, y_dirac=s)
        pp=ppde.getSolution()
        uu=loc(pp)

        # arguments for DcRes
        current = 10.
        sourceInfo = [ "sl0",  "sl1" ]
        sampleTags = [ ("sr0", "sr1") ]

        sigmaPrimary=7.
        phiPrimary=pp*current*sigma0/sigmaPrimary

        uuscale=1-current*sigma0/sigmaPrimary
        delphi_in = [ (uu[1]-uu[0]) * uuscale]

        acw=DcRes(domain, loc, delphi_in, sampleTags,  phiPrimary, sigmaPrimary)

        self.assertLess(Lsup(phiPrimary-acw.getPrimaryPotential()), 1.e-10 * Lsup(acw.getPrimaryPotential()))

        SIGMA=10. # matches current
        args0=acw.getArguments(SIGMA)
        p=args0[0]
        u=args0[1]

        # true secondary potential
        pps=pp-phiPrimary
        self.assertLess(Lsup(p-pps), 1.e-6*Lsup(pps))

        # test return values at electrodes:
        self.assertLess(abs(u[0]-uu[0]*uuscale), 1.e-6 * abs(uu[0]*uuscale))
        self.assertLess(abs(u[1]-uu[1]*uuscale), 1.e-6 * abs(uu[1]*uuscale))

        # this sould be zero
        dd=acw.getDefect(SIGMA, *args0)
        self.assertTrue( dd >= 0.)
        self.assertTrue( dd <= 1e-7 )

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
   def test_StrongCoupled4(self):
        x=Solution(self.domain).getX()
        # --- set exact solution ----
        u_ex=Data(1.,(4,),Solution(self.domain))
        g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
        A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
        for i in range(self.domain.getDim()):
           u_ex[0]+= 1*(i+1)*x[i]
           g_ex[0,i]=1*(i+1)
           u_ex[1]+= 2*(i+1)*x[i]
           g_ex[1,i]=2*(i+1)
           u_ex[2]+= 3*(i+1)*x[i]
           g_ex[2,i]=3*(i+1)
           u_ex[3]+= 4*(i+1)*x[i]
           g_ex[3,i]=4*(i+1)
           A[0,i,0,i]=1
           A[1,i,1,i]=1
           A[2,i,2,i]=1
           A[3,i,3,i]=1

        a=-1./4.*0.99

        Y=Data(0.,(4,),Function(self.domain))
        Y[0]=     u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
        Y[1]= a * u_ex[0]+     u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
        Y[2]= a * u_ex[0]+ a * u_ex[1]+     u_ex[2]+ a * u_ex[3]
        Y[3]= a * u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+     u_ex[3]
        # create PDE:
        pde=LinearPDE(self.domain,numEquations=4)
        pde.setSymmetryOn()
        mask=whereZero(x[0])*[1,1,1,1]
        pde.setValue(r=u_ex,q=mask)
        pde.setValue(A=A,
                     D=kronecker(4)+a*(numpy.ones((4,4))-kronecker(4)),
                     Y=Y,
                     y=matrixmult(g_ex,self.domain.getNormal()))
        # -------- get the solution ---------------------------
        pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
        pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
        pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
        pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
        if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
        if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
        if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
        u=pde.getSolution()
        # -------- test the solution ---------------------------
        error=Lsup(u-u_ex)/Lsup(u_ex)
        self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
    def getPotentialNumeric(self):
        """
        Returns 3 list each made up of a number of list containing primary, secondary and total
        potentials diferences. Each of the lists contain a list for each value of n.
        """
        primCon=self.primaryConductivity
        coords=self.domain.getX()
        tol=1e-8
        pde=LinearPDE(self.domain, numEquations=1)
        pde.getSolverOptions().setTolerance(tol)
        pde.setSymmetryOn()
        primaryPde=LinearPDE(self.domain, numEquations=1)
        primaryPde.getSolverOptions().setTolerance(tol)
        primaryPde.setSymmetryOn()
        DIM=self.domain.getDim()
        x=self.domain.getX()
        q=whereZero(x[DIM-1]-inf(x[DIM-1]))
        for i in xrange(DIM-1):
            xi=x[i]
            q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
        A = self.secondaryConductivity * kronecker(self.domain)
        APrimary = self.primaryConductivity * kronecker(self.domain)
        pde.setValue(A=A,q=q)
        primaryPde.setValue(A=APrimary,q=q)

        delPhiSecondaryList = []
        delPhiPrimaryList = []
        delPhiTotalList = []
        for i in range(1,self.n+1): # 1 to n
            maxR = self.numElectrodes - 1 - (2*i) #max amount of readings that will fit in the survey
            delPhiSecondary = []
            delPhiPrimary = []
            delPhiTotal = []
            for j in range(maxR):
                y_dirac=Scalar(0,DiracDeltaFunctions(self.domain))
                y_dirac.setTaggedValue(self.electrodeTags[j],self.current)
                y_dirac.setTaggedValue(self.electrodeTags[j + ((2*i) + 1)],-self.current)
                self.sources.append([self.electrodeTags[j], self.electrodeTags[j + ((2*i) + 1)]])
                primaryPde.setValue(y_dirac=y_dirac)
                numericPrimaryPot = primaryPde.getSolution()
                X=(primCon-self.secondaryConductivity) * grad(numericPrimaryPot)
                pde.setValue(X=X)
                u=pde.getSolution()
                loc=Locator(self.domain,[self.electrodes[j+i],self.electrodes[j+i+1]])
                self.samples.append([self.electrodeTags[j+i],self.electrodeTags[j+i+1]])
                valPrimary=loc.getValue(numericPrimaryPot)
                valSecondary=loc.getValue(u)
                delPhiPrimary.append(valPrimary[1]-valPrimary[0])
                delPhiSecondary.append(valSecondary[1]-valSecondary[0])
                delPhiTotal.append(delPhiPrimary[j]+delPhiSecondary[j])
            delPhiPrimaryList.append(delPhiPrimary)
            delPhiSecondaryList.append(delPhiSecondary)
            delPhiTotalList.append(delPhiTotal)

        self.delPhiPrimaryList=delPhiPrimaryList
        self.delPhiSecondaryList=delPhiSecondaryList
        self.delPhiTotalList = delPhiTotalList
        return [delPhiPrimaryList, delPhiSecondaryList, delPhiTotalList]

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
def runTest(dom, n=1, a=0, b=0):
   print("================== TEST : n= %s a=%s b=%s ================"%(n,a,b))
   DIM=dom.getDim()
   normal=dom.getNormal()
   mypde=LinearPDE(dom, numEquations=n, numSolutions=n)
   x=dom.getX()
   A=mypde.createCoefficient("A")
   D=mypde.createCoefficient("D")
   Y=mypde.createCoefficient("Y")
   y=mypde.createCoefficient("y")
   q=mypde.createCoefficient("q")
   if n==1:
      u_ref=Scalar(0.,Solution(dom))
      for j in range(DIM):
         q+=whereZero(x[j])
         A[j,j]=1
      y+=sin(sqrt(2.))*normal[0]
      Y+=b*x[0]*sin(sqrt(2.))
      D+=b
      u_ref=x[0]*sin(sqrt(2.))
   else:
      u_ref=Data(0.,(n,),Solution(dom))
      for i in range(n):
          for j in range(DIM):
             q[i]+=whereZero(x[j])
             A[i,j,i,j]=1
             if j == i%DIM: y[i]+=sin(i+sqrt(2.))*normal[j]
          for j in range(n):
                if i==j:
                   D[i,i]=b+(n-1)*a
                   Y[i]+=b*x[i%DIM]*sin(i+sqrt(2.))
                else:
                   D[i,j]=-a
                   Y[i]+=a*(x[i%DIM]*sin(i+sqrt(2.))-x[j%DIM]*sin(j+sqrt(2.)))
          u_ref[i]=x[i%DIM]*sin(i+sqrt(2.))
          
#    - u_{j,ii} + b*u_i+ a*sum_{k<>j}  (u_i-u_k) = F_j
          
   mypde.setValue(A=A, D=D, y=y, q=q, r=u_ref, Y=Y)
   mypde.getSolverOptions().setVerbosityOn()
   mypde.getSolverOptions().setTolerance(TOL)
   mypde.setSymmetryOn()
   u=mypde.getSolution()
   error=Lsup(u-u_ref)/Lsup(u_ref)
   print("error = ",error)
   if error > 10*TOL: print("XXXXXXXXXX Error to large ! XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX")

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
    def gzpot(p, y, x, *args):
        #rho, rhox, rhoy, R = p
        rhox=args[0]/2.; rhoy=args[1]/2.
        rho, R, z =p
        #Domain related.
        mx = args[0]; my = args[1];
        #PDE related
        G=6.67300*10E-11

        #DOMAIN CONSTRUCTION
        domain=ReadGmsh('data/example10m/example10m.msh',2)
        domx=Solution(domain).getX()
        mask=wherePositive(R-length(domx-rholoc))
        rhoe=rho*mask
        kro=kronecker(domain)

        q=whereZero(domx[1]-my)+whereZero(domx[1])+whereZero(domx[0])+whereZero(domx[0]-mx)
        #ESCRIPT PDE CONSTRUCTION
        mypde=LinearPDE(domain)
        mypde.setValue(A=kro,Y=4.*np.pi*G*rhoe,q=q,r=0.0)
        mypde.setSymmetryOn()
        sol=mypde.getSolution()

        g_field=grad(sol) #The graviational accelleration g.
        g_fieldz=g_field*[0,1] #The vertical component of the g field.
        gz=length(g_fieldz) #The magnitude of the vertical component.

        #MODEL SIZE SAMPLING
        sol_escgz=[]
        sol_escx=[]
        for i in range(0,len(x)):
            sol_escgz.append([x[i],rhoy+z])

        sample=[] # array to hold values
        rec=Locator(gz.getFunctionSpace(),sol_escgz) #location to record
        psol=rec.getValue(gz)

        err = np.sum((np.array(y) - np.array(psol))**2.)
        print("Lsup= ",Lsup(np.array(psol)-np.array(sol_angz))/Lsup(np.array(psol)))
        return err

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
class StokesProblem(SaddlePointProblem):
      """
      simple example of saddle point problem
      """
      def __init__(self,domain,debug=False):
         super(StokesProblem, self).__init__(self,debug)
         self.domain=domain
         self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
         self.__pde_u.setSymmetryOn()

         self.__pde_p=LinearPDE(domain)
         self.__pde_p.setReducedOrderOn()
         self.__pde_p.setSymmetryOn()

      def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1):
         self.eta=eta
         A =self.__pde_u.createCoefficientOfGeneralPDE("A")
         for i in range(self.domain.getDim()):
           for j in range(self.domain.getDim()):
             A[i,j,j,i] += self.eta
             A[i,j,i,j] += self.eta
         self.__pde_p.setValue(D=1/self.eta)
         self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f)

      def inner(self,p0,p1):
         return integrate(p0*p1,Function(self.__pde_p.getDomain()))

      def solve_f(self,u,p,tol=1.e-8):
         self.__pde_u.setTolerance(tol)
         g=grad(u)
         self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain()))
         return  self.__pde_u.getSolution()

      def solve_g(self,u,tol=1.e-8):
         self.__pde_p.setTolerance(tol)
         self.__pde_p.setValue(X=-u) 
         dp=self.__pde_p.getSolution()
         return  dp

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
class SlippingFault(SaddlePointProblem):
      """
      simple example of saddle point problem
      """
      def __init__(self,domain):
         super(SlippingFault, self).__init__(self)
         self.domain=domain
         self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
         self.__pde_u.setSymmetryOn()

      def initialize(self,density=1.,lmbd=1., mu=1., traction=Data(),fixed_u_mask=Data(), slip=0.):
         d=self.domain.getDim()
         self.slip=slip
         A =self.__pde_u.createCoefficientOfGeneralPDE("A")
         for i in range(self.domain.getDim()):
           for j in range(self.domain.getDim()):
             A[i,j,j,i] += mu
             A[i,j,i,j] += mu
             A[i,i,j,j] += lmbd
         self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=-kronecker(Function(self.domain))[d-1]*g*density,y=traction)

      def inner(self,p0,p1):
         return integrate(inner(p0,p1),FunctionOnContactZero(self.domain))

      def solve_f(self,u,p,tol=1.e-8):
         self.__pde_u.setTolerance(tol)
         self.__pde_u.setValue(y_contact=-p)
         # print "p:",inf(p),sup(p)
         # print "u:",inf(u),sup(u)
         self.__pde_u.setValue(y_contact=-p)
         return  self.__pde_u.getSolution()

      def solve_g(self,u,tol=1.e-8):
         dp=Vector(0.,FunctionOnContactZero(self.domain))
         h=FunctionOnContactZero(self.domain).getSize()
         # print jump(u)-self.slip
         dp[0]=(self.slip[0]-jump(u[0]))*lam_mu/h
         dp[1]=(self.slip[1]-jump(u[1]))*lam_mu/h
         dp[2]=(self.slip[2]-jump(u[2]))*lam_mu/h
         return  dp

# 需要导入模块: from esys.escript.linearPDEs import LinearPDE [as 别名]
# 或者: from esys.escript.linearPDEs.LinearPDE import setValue [as 别名]
 # material parameter
 k=1; rhocp=100; 
 # bottom temperature:
 T_bot=100
 # generate domain:
 mydomain=Rectangle(l0=L0,l1=L1,n0=20,n1=20)
 x=mydomain.getX()
 # set boundray temperature:
 T_D=T_bot/L1*(L1-x[1])
 # set heat source:
 Q=Qc*whereNegative(length(x-xc)-r)
 # generate domain:
 mypde=LinearPDE(mydomain)
 mypde.setSymmetryOn()
 # set PDE coefficients:
 mypde.setValue(A=dt*k*kronecker(mydomain), D=dt*rhocp, 
                 r=T_D, q=whereZero(x[1])+whereZero(x[1]-L1))
 # initial temperature
 T=T_D 
 # step counter and time marker:
 N=0; t=0
 # stop when t_end is reached:
 while t<t_end:
     print("time step %d, t=%s"%(N,t))
     # update PDE coefficient:
     mypde.setValue(Y=dt*rhocp*T+dt*Q)
     # new temperature:
     T=mypde.getSolution()
     # save as VTK for visualisation:
     saveVTK("u.%s.vtu"%N,T=T)
     # increase counter and marker:
     N+=1; t+=dt