Module regpy.operators.ngsolve
PDE forward operators using NGSolve
Classes
class NGSolveOperator (domain, codomain, linear=False)-
Expand source code
class NGSolveOperator(Operator): def __init__(self, domain, codomain, linear = False): super().__init__(domain = domain, codomain = codomain, linear = linear) self.gfu_read_in = ngs.GridFunction(self.domain.fes) '''Reads in a coefficient vector of the domain and interpolates in the codomain. The result is saved in gfu''' def _read_in(self, vector, gfu,definedonelements=None): """Read in of a numpy array into a ngsolve grid function. Note You can also read into ngsolve LinearForm. Parameters ---------- vector : numpy.array Numpy array to be put into the `GridFunction` gfu : ngsolve.GridFunction or ngsolve.LinearForm ngsolve element to be written into. definedonelements : pyngcore.pyngcore.BitArray, optional BitArray representing the finite elements to be projected onto, by default None """ gfu.vec.FV().NumPy()[:] = vector if definedonelements is not None: ngs.Projector(definedonelements, range=True).Project(gfu.vec) def _solve_dirichlet_problem(self, bf, lf, gf, prec): r"""Solves the problem \begin{align*} b(u,v) = f(v) \;\forall v\; test\; functions,\\ u|_\Gamma = g \end{align*} The boundary values are given by what `gf` has as values on the boundary. Parameters ---------- bf : ngs.BilinearForm the bilinear form of the problem. lf : ngs.LinearFrom the linear form of the problem. gf : ngs.GridFunction The grid functions on which to solve the solution will be put into these and they have to satisfy the boundary condition that you want. prec : ngs.Preconditioner preconditioner to be used with ngsolve. """ gf.vec.data += bf.mat.Inverse(freedofs=self.codomain.fes.FreeDofs()) * (lf.vec - bf.mat * gf.vec)Ancestors
Subclasses
- EIT
- ElasticAnisotrop
- MapCubic
- MapIsotropic
- MapTransverse
- ProjectToBoundary
- ReactionNeumann
- SecondOrderEllipticCoefficientPDE
Inherited members
class SecondOrderEllipticCoefficientPDE (domain, sol_domain, bdr_val=None, a_bdr_val=None, adj_lf=True, cu_lf=True)-
Expand source code
class SecondOrderEllipticCoefficientPDE(NGSolveOperator): r"""Provides a general setup for the forward problems mapping PDE coefficients to their solutions. That is we assume that as variational formulation one can use \[ \forall v: b_0(u,v) + b_a(u,v) = F(v) \] with \(b_0\) a bilinear form independent of the coefficient \(a\) and \(F\) some linear form. Furthermore \(b_a\) the bilinear form that depends on \(a\) has to be linear in \(a\) so that one can define a bilinear form \(c_u(a,v)=b_a(u,v)\). More over we may assume some Dirichlet boundary conditions on \(u\).That is \[ F: a \mapsto u \] The Frechét derivative \(F'[a]h\) in direction \(h\) is then given as the variational solution \(u'\) to \[ \forall v: b_0(u',v) + b_a(u',v) = -c_u(h,v) \] with \(u=F(a)\). That is \[ F'[a]: h \mapsto u' \] It's adjoint \(F'[a]^\ast g\) is given as the linear form \(-c_u(\cdot,w)\) for \(u=F(a)\) and \(w\) solving the problem \[ \forall v: b_0(v,w) + b_a(v,w) = <g,v>. \] That is \[ F'[a]^\ast: g \mapsto -c_u(\cdot,w) \] Notes ----- A subclass implemented by a user has to at least implement the subroutine `_bf` which is the implementation of the bilinear form depending on the coefficient \(a\). Optional can the independent bilinear form \(b_0\) be implemented as formal integrator in `_bf_0` and the linear form can be implemented in `_lf`. Parameters ---------- domain : NgsSpace The NgsSpace on which the coefficients defined are. sol_domain : NgsSpace The NgsSpace on which the PDE solutions defined are. bdr_val : array type, optional Boundary value of the PDE solution of the forward evaluation, by default None a_bdr_val : array type, optional Boundary value of the coefficients, by default None """ def __init__(self, domain, sol_domain, bdr_val = None, a_bdr_val = None,adj_lf=True,cu_lf=True): super().__init__(domain, sol_domain, linear = False) self.gfu_a=ngs.GridFunction(self.domain.fes) self.gfu_h=ngs.GridFunction(self.domain.fes) self.gfu_a_bdr = ngs.GridFunction(self.domain.fes) if a_bdr_val is not None: assert self.domain.bdr is not None assert self.domain.is_on_boundary(a_bdr_val) self._read_in(a_bdr_val,self.gfu_a_bdr) self.gfu_deriv=ngs.GridFunction(self.codomain.fes) self.gfu_adj_help=ngs.GridFunction(self.codomain.fes) if bdr_val is not None and bdr_val in self.codomain: self.gfu_eval=self.codomain.to_ngs(bdr_val) else: self.gfu_eval=ngs.GridFunction(self.codomain.fes) self.adj_lf=adj_lf self.cu_lf=cu_lf self.u_a, self.v_a = self.domain.fes.TnT() self.u, self.v = self.codomain.fes.TnT() def _eval(self, a, differentiate=False, adjoint_derivative=False): self._read_in(a, self.gfu_a,definedonelements=self.domain.fes.FreeDofs()) #self.gfu_a.vec.data = self.gfu_a.vec + self.gfu_a_bdr.vec self.bf_mat = ngs.BilinearForm(self.codomain.fes) self.bf_mat += self._bf(self.gfu_a,self.u,self.v) if self._bf_0() is not None: self.bf_mat += self._bf_0() self.bf_mat.Assemble() self.bf_mat_inv = self.bf_mat.mat.Inverse(freedofs=self.codomain.fes.FreeDofs()) self.gfu_eval.vec.data += self.bf_mat_inv * (self._lf().vec - self.bf_mat.mat * self.gfu_eval.vec) return self.gfu_eval.vec.FV().NumPy().copy() def _derivative(self, h): self._read_in(h, self.gfu_h,definedonelements=self.domain.fes.FreeDofs()) lf = ngs.LinearForm(self.codomain.fes) lf = self._c_u(h) self.gfu_deriv.vec.data += self.bf_mat_inv * (-lf.vec - self.bf_mat.mat * self.gfu_deriv.vec) return self.gfu_deriv.vec.FV().NumPy().copy() def _adjoint(self, g): lf = ngs.LinearForm(self.codomain.fes).Assemble() self._read_in(g,lf) self.gfu_adj_help.vec.data += self.bf_mat.mat.CreateTranspose().Inverse(freedofs=self.codomain.fes.FreeDofs()) * (lf.vec - self.bf_mat.mat.CreateTranspose() * self.gfu_adj_help.vec) if self.adj_lf: lf_adj = ngs.LinearForm(self.domain.fes) lf_adj += -1*self._bf(self.v_a,self.gfu_eval,self.gfu_adj_help) lf_adj.Assemble() ngs.Projector(self.domain.fes.FreeDofs(), range=True).Project(lf_adj.vec) return lf_adj.vec.FV().NumPy().copy() else: grid_adj=ngs.GridFunction(self.domain.fes) strain1= ngs.Sym(ngs.Grad(self.gfu_eval)) strain2=ngs.Sym(ngs.Grad(self.gfu_adj_help)) grid_adj=self._eval_adj(strain1,strain2) ngs.Projector(self.domain.fes.FreeDofs(), range=True).Project(grid_adj.vec) return grid_adj.vec.FV().NumPy().copy() def _bf(self,a,u,v): r"""Implementation of \(b_a\) as `ngsolve.comp.SumOfIntegrals` that is something similar to `a*ngs.grad(u)*ngs.grad(v)*ngs.dx` where `u` ist used as trial functions and `v` as test functions. This method has to be implemented by Parameters ---------- a : ngsolve.CoefficientFunction or ngsolve.GridFunction the coefficient in the PDE u : ngsolve.comp.ProxyFunction Trial functions for PDE v : ngsolve.comp.ProxyFunction Test functions for PDE Returns ------ ngsolve.comp.SumOfIntegrals the formal Integration formula of the bilinear form depending on the coefficient """ raise NotImplementedError def _bf_0(self): r"""Implementation of \(b_0\) as `ngsolve.comp.SumOfIntegrals` is an optional method to be overwritten with subclasses. Returns ------- ngsolve.comp.SumOfIntegrals the formal Integration formula of the bilinear form independent of the coefficient """ return None def _lf(self): r"""The Linear form of the PDE \(F\) implemented as a fixed Linear form. Note that the Linear form has to be defined on the `codomain` as this is the domain of the solution of the PDE. By default this is the empty Linear form. Returns ------ ngsolve.Linearform Linear form of the PDE \(F\) as `ngsolve.LinearForm` """ return ngs.LinearForm(self.codomain.fes).Assemble() def _c_u(self,h,c_u_help=False): self._read_in(h, self.gfu_h,definedonelements=self.domain.fes.FreeDofs()) lf = ngs.LinearForm(self.codomain.fes) lf += self._bf(self.gfu_h,self.gfu_eval,self.v) return lf.Assemble() def _eval_adj(self,u,v): return NoneProvides a general setup for the forward problems mapping PDE coefficients to their solutions. That is we assume that as variational formulation one can use \forall v: b_0(u,v) + b_a(u,v) = F(v) with b_0 a bilinear form independent of the coefficient a and F some linear form.
Furthermore b_a the bilinear form that depends on a has to be linear in a so that one can define a bilinear form c_u(a,v)=b_a(u,v). More over we may assume some Dirichlet boundary conditions on u.That is F: a \mapsto uThe Frechét derivative F'[a]h in direction h is then given as the variational solution u' to \forall v: b_0(u',v) + b_a(u',v) = -c_u(h,v) with u=F(a). That is F'[a]: h \mapsto u'
It's adjoint F'[a]^\ast g is given as the linear form -c_u(\cdot,w) for u=F(a) and w solving the problem \forall v: b_0(v,w) + b_a(v,w) = <g,v>. That is F'[a]^\ast: g \mapsto -c_u(\cdot,w)
Notes
A subclass implemented by a user has to at least implement the subroutine
_bfwhich is the implementation of the bilinear form depending on the coefficient a. Optional can the independent bilinear form b_0 be implemented as formal integrator in_bf_0and the linear form can be implemented in_lf.Parameters
domain:NgsSpace- The NgsSpace on which the coefficients defined are.
sol_domain:NgsSpace- The NgsSpace on which the PDE solutions defined are.
bdr_val:array type, optional- Boundary value of the PDE solution of the forward evaluation, by default None
a_bdr_val:array type, optional- Boundary value of the coefficients, by default None
Ancestors
Inherited members
class ElasticAnisotrop (domain, sol_domain, rhs, bdr_val=None, bdr_N=None, bdr_N_val=None)-
Expand source code
class ElasticAnisotrop(NGSolveOperator): r'''Operator for solving the general elasticity problem with 4th order material tensor tensor is symmetrized, but corresponding FE-space has dimension d x d x d x d. ''' def __init__(self, domain, sol_domain, rhs, bdr_val = None, bdr_N=None, bdr_N_val=None): super().__init__(domain, sol_domain, linear = False) #gridfunctions in domain of material tensor #gridfunction for material tensor self.gfu_a=ngs.GridFunction(self.domain.fes) #gridfunction for h-value in derivative self.gfu_h=ngs.GridFunction(self.domain.fes) #gridfunction for adjoint self.gfu_adj=ngs.GridFunction(self.domain.fes) #dimension of the finite element mesh self.dim=self.codomain.fes.mesh.dim #gridfunctions in codomain of finite eleament solutions #derivative self.gfu_deriv=ngs.GridFunction(self.codomain.fes) #solution of adjoint FE problem for adjoint of derivative self.gfu_adj_help=ngs.GridFunction(self.codomain.fes) #rhs of adjoint FE-problem self.gfu_adj_rhs=ngs.GridFunction(self.codomain.fes) #Neumann boundary conditions self.gfu_neumann=ngs.GridFunction(self.codomain.fes) #solution of forward problem self.gfu_eval=ngs.GridFunction(self.codomain.fes) #Dirichlet boundary conditions if bdr_val is not None and bdr_val in self.codomain: self.gfu_eval.Set(bdr_val) else: if self.dim==2: self.gfu_eval.Set((0,0)) elif self.dim==3: self.gfu_eval.Set((0,0,0)) else: raise NotImplementedError #Neumann boundary conditions self.bdr_N=bdr_N if self.bdr_N is not None and bdr_N_val is not None: #assert bdr_N_val in self.codomain self.gfu_neumann.Set(bdr_N_val) #right hand side of forward problem self.rhs=rhs #Dirichlet boundary conditions for derivative and adjoint if self.dim==2: self.gfu_deriv.Set((0,0)) self.gfu_adj_help.Set((0,0)) elif self.dim==3: self.gfu_deriv.Set((0,0,0)) self.gfu_adj_help.Set((0,0,0)) else: raise NotImplementedError #test and trial functions self.u, self.v = self.codomain.fes.TnT() def _eval(self, a, differentiate=False, adjoint_derivative=False): #a: fourth order tensor in terms of flat numpy array self._read_in(a, self.gfu_a,definedonelements=self.domain.fes.FreeDofs()) self.bf_mat = ngs.BilinearForm(self.codomain.fes) self.prec=ngs.Preconditioner(self.bf_mat, "local") #use matrix-vector multiplication for tensor contraction dim2=self.dim*self.dim a_mat=ngs.CF(self.gfu_a,dims=(dim2,dim2)) u_vec=ngs.CF(ngs.Sym(ngs.Grad(self.u)),dims=(dim2,)) v_vec=ngs.CF(ngs.Sym(ngs.Grad(self.v)),dims=(dim2,)) sigu=a_mat*u_vec sigv=a_mat*v_vec #Bilinear form to be used in other methods self.bf_mat += (ngs.InnerProduct(sigu,ngs.Sym(ngs.Grad(self.v)))+ngs.InnerProduct(sigv,ngs.Sym(ngs.Grad(self.u))))*0.5 * ngs.dx with ngs.TaskManager(): self.bf_mat.Assemble() #Linear form lf = ngs.LinearForm(self.codomain.fes) lf += ngs.InnerProduct(self.rhs, self.v)* ngs.dx if self.bdr_N is not None: lf += ngs. InnerProduct(self.gfu_neumann,self.v)* ngs.ds(self.bdr_N) with ngs.TaskManager(): lf.Assemble() #Solution of bdr problem, inverse of bilinear form will be used in other methods self.bf_mat_inv = self.bf_mat.mat.Inverse(freedofs=self.codomain.fes.FreeDofs()) self.gfu_eval.vec.data += self.bf_mat_inv * (lf.vec - self.bf_mat.mat * self.gfu_eval.vec) #result of application of forward operator return self.gfu_eval.vec.FV().NumPy().copy() def _derivative(self, h): #h: fourth order tensor in terms of flat numpy array self._read_in(h, self.gfu_h,definedonelements=self.domain.fes.FreeDofs()) lf = ngs.LinearForm(self.codomain.fes) #use matrix-vector multiplication for tensor contraction dim2=self.dim*self.dim a_mat=ngs.CF(self.gfu_h,dims=(dim2,dim2)) u_vec=ngs.CF(ngs.Sym(ngs.Grad(self.gfu_eval)),dims=(dim2,)) v_vec=ngs.CF(ngs.Sym(ngs.Grad(self.v)),dims=(dim2,)) sigu=a_mat*u_vec sigv=a_mat*v_vec #rhs for pde determining Frechet derivative lf +=(ngs.InnerProduct(sigu,ngs.Sym(ngs.Grad(self.v)))+ngs.InnerProduct(sigv,ngs.Sym(ngs.Grad(self.gfu_eval))))*0.5 * ngs.dx lf.Assemble() #solution of pde for frechet derivative self.gfu_deriv.vec.data += self.bf_mat_inv * (-lf.vec - self.bf_mat.mat * self.gfu_deriv.vec) return self.gfu_deriv.vec.FV().NumPy().copy() def _adjoint(self, g): #solve adjoint forward PDE lf = ngs.LinearForm(self.codomain.fes) self._read_in(g,self.gfu_adj_rhs) lf +=ngs.InnerProduct(self.gfu_adj_rhs, self.v) * ngs.dx lf.Assemble() self.gfu_adj_help.vec.data += self.bf_mat.mat.CreateTranspose().Inverse(freedofs=self.codomain.fes.FreeDofs()) * (lf.vec - self.bf_mat.mat.CreateTranspose() * self.gfu_adj_help.vec) gfu_adj=ngs.GridFunction(self.domain.fes) #adjoint is symmetrized dyadic product strain1= ngs.Sym(ngs.Grad(self.gfu_eval)) strain2=ngs.Sym(ngs.Grad(self.gfu_adj_help)) list_adj=[] dim2=self.dim*self.dim for i in range(dim2): for j in range (dim2): sijkl=-0.5*(strain1[i]*strain2[j]+strain2[i]*strain1[j]) list_adj.append(sijkl) #direct transfer to gfu_adj does not work because sijkl are gridfunctions help1=ngs.CF(tuple(list_adj)) gfu_adj.Set(help1) ngs.Projector(self.domain.fes.FreeDofs(), range=True).Project(gfu_adj.vec) return gfu_adj.vec.FV().NumPy().copy()Operator for solving the general elasticity problem with 4th order material tensor tensor is symmetrized, but corresponding FE-space has dimension d x d x d x d.
Ancestors
Inherited members
class MapIsotropic (domain, codomain, nu=0.3)-
Expand source code
class MapIsotropic(NGSolveOperator): r"""maps physical input parameters to subspace in material tensor space isotropic elasticity """ def __init__(self, domain, codomain,nu=0.3 ): super().__init__(domain = domain, codomain = codomain, linear = True) #material parameters self.gfu_a=ngs.GridFunction(self.domain.fes) #4th order material tensor self.gfu_tens=ngs.GridFunction(self.codomain.fes) #container for adjoint self.gfu_map=ngs.GridFunction(self.codomain.fes) #number of free parameters self.dim_domain=1 #material tensor has 16 components in 2D case and 81 in 3D self.dim=self.codomain.fes.mesh.dim self.dim_codomain=self.dim**4 self.nu=nu #container for material data self.list_mat=[] for i in range(self.dim_codomain): self.list_mat.append(ngs.CF(0)) #assign generalized Lame constants to material tensor if self.dim_codomain==16: self.list_lam1=[0,15] self.list_lam4=[3,12] self.list_mu=[5,6,9,10] else: self.list_lam1=[0,40,80] self.list_lam4=[4,8,36,44,72,76] self.list_mu=[10,12,20,24,28,30,50,52,56,60,68,70] def _eval(self,argument): self._read_in(argument, self.gfu_a) #elastic constants E=self.gfu_a lam1=E*(1-self.nu)/((1+self.nu)*(1-2*self.nu)) lam4=E*self.nu/((1+self.nu)*(1-2*self.nu)) mu=E/(2*(1+self.nu)) #definition of material tensor for i in self.list_lam1: self.list_mat[i]=lam1 for i in self.list_lam4: self.list_mat[i]=lam4 for i in self.list_mu: self.list_mat[i]=mu #convert to gridfunction for further use self.gfu_tens.Set(tuple(self.list_mat)) #flat numpy-array from gridfunction return self.gfu_tens.vec.FV().NumPy().copy() def _adjoint(self,argument): self._read_in(argument, self.gfu_map) #material tensor/E gfu_help=ngs.GridFunction(self.codomain.fes) help1=self._eval(self.domain.ones()) self._read_in(help1, gfu_help) #adjoint help2=ngs.InnerProduct(self.gfu_map,gfu_help) self.gfu_a.Set(help2) return self.gfu_a.vec.FV().NumPy().copy()maps physical input parameters to subspace in material tensor space isotropic elasticity
Ancestors
Inherited members
class MapCubic (domain, codomain, nu=0.3)-
Expand source code
class MapCubic(NGSolveOperator): r"""projects physical input parameters to subspace in material tensor space cubic elasticity """ def __init__(self, domain, codomain,nu=0.3 ): super().__init__(domain = domain, codomain = codomain, linear = True) #grid function for physical parameters self.gfu_a=ngs.GridFunction(self.domain.fes) #grid function for material tensor self.gfu_tens=ngs.GridFunction(self.codomain.fes) #grid function for input of adjoint self.gfu_map=ngs.GridFunction(self.codomain.fes) self.nu=nu #material tensor has 16 components in 2D case and 81 in 3D self.dim=self.codomain.fes.mesh.dim self.dim_codomain=self.dim**4 #number of free parameters self.dim_domain = 2 #container for material data self.list_mat=[] for i in range(self.dim_codomain): self.list_mat.append(ngs.CF(0)) #assign generalized Lame constants to material tensor if self.dim_codomain==16: self.list_lam1=[0,15] self.list_lam4=[3,12] self.list_mu=[5,6,9,10] else: self.list_lam1=[0,40,80] self.list_lam4=[4,8,36,44,72,76] self.list_mu=[10,12,20,24,28,30,50,52,56,60,68,70] #parameter specitic components of material tensor help1=self.codomain.to_ngs(self._eval(self.domain.ones())) self.help_E,self.help_mu=self._split_map(help1) def _split_map(self,help1): #split projector map_E=[] map_mu=[] for i in range(self.dim_codomain): if i in self.list_mu: map_mu.append(help1.components[i]) map_E.append(ngs.CF(0)) elif i in self.list_lam1 or i in self.list_lam4: map_E.append(help1.components[i]) map_mu.append(ngs.CF(0)) else: map_E.append(ngs.CF(0)) map_mu.append(ngs.CF(0)) help_E=ngs.GridFunction(self.codomain.fes) help_mu=ngs.GridFunction(self.codomain.fes) help_E.Set(tuple(map_E)) help_mu.Set(tuple(map_mu)) return help_E, help_mu def _eval(self,argument): self._read_in(argument, self.gfu_a) #elastic constants self._read_in(argument, self.gfu_a) #elastic constants E=self.gfu_a.components[0] mu=self.gfu_a.components[1] lam1=E*(1-self.nu)/((1+self.nu)*(1-2*self.nu)) lam4=E*self.nu/((1+self.nu)*(1-2*self.nu)) #mu=E/(2*(1+self.nu)) #definition of material tensor for i in self.list_lam1: self.list_mat[i]=lam1 for i in self.list_lam4: self.list_mat[i]=lam4 for i in self.list_mu: self.list_mat[i]=mu #convert to gridfunction for further use self.gfu_tens.Set(tuple(self.list_mat)) #flat numpy-array from gridfunction return self.gfu_tens.vec.FV().NumPy().copy() def _adjoint(self,argument): self._read_in(argument, self.gfu_map) adj_E=ngs.InnerProduct(self.gfu_map,self.help_E) adj_mu=ngs.InnerProduct(self.gfu_map, self.help_mu) help2=(adj_E,adj_mu) self.gfu_a.Set(help2) return self.gfu_a.vec.FV().NumPy().copy()projects physical input parameters to subspace in material tensor space cubic elasticity
Ancestors
Inherited members
class MapTransverse (domain, codomain, nu_A=0.2, nu_T=0.3, scale=0.3)-
Expand source code
class MapTransverse(NGSolveOperator): r"""maps physical input parameters to subspace in material tensor space transverse isotropic elastice with isotropic plane 23 """ def __init__(self, domain, codomain,nu_A=0.2,nu_T=0.3,scale=0.3): super().__init__(domain = domain, codomain = codomain, linear = True) self.gfu_a=ngs.GridFunction(self.domain.fes) self.gfu_tens=ngs.GridFunction(self.codomain.fes) self.gfu_map=ngs.GridFunction(self.codomain.fes) self.nu_T=nu_T self.nu_A=nu_A self.scale=scale #material tensor self.dim=self.codomain.fes.mesh.dim self.dim_codomain=self.dim**4 #number of free parameters self.dim_domain = len(self.domain.fes.components) #contaimer for material data self.list_mat=[] for i in range(self.dim_codomain): self.list_mat.append(ngs.CF(0)) #assigns generalized Lame constants to material tensor if self.dim_codomain==16: self.list_lam4=[3,12] self.list_mu2=[5,6,9,10] #need dummy lists for adjoint self.list_lam6=[] self.list_lam2=[15] if self.dim_codomain==81: self.list_lam2=[40,80] self.list_lam4=[4,8,36,72] self.list_lam6=[44,76] self.list_mu1=[50,52,68,70] self.list_mu2=[10,12,20,24,28,30,56,60] help1=self.codomain.to_ngs(self._eval(self.domain.ones())) self.help_EA,self.help_ET,self.help_mu2=self._split_map(help1) def _split_map(self,help1): #split material operator map_EA=[help1.components[0]] map_ET=[0] map_mu2=[0] for i in range(1,self.dim_codomain): if i in self.list_mu2: map_mu2.append(help1.components[i]) map_ET.append(ngs.CF(0)) map_EA.append(ngs.CF(0)) elif self.dim==2 and (i in self.list_lam2 or i in self.list_lam4 or i in self.list_lam6): map_mu2.append(ngs.CF(0)) map_ET.append(help1.components[i]) map_EA.append(ngs.CF(0)) elif self.dim==3 and (i in self.list_lam2 or i in self.list_lam4 or i in self.list_lam6 or i in self.list_mu1): map_mu2.append(ngs.CF(0)) map_ET.append(help1.components[i]) map_EA.append(ngs.CF(0)) else: map_mu2.append(ngs.CF(0)) map_ET.append(ngs.CF(0)) map_EA.append(ngs.CF(0)) help_ET=ngs.GridFunction(self.codomain.fes) help_EA=ngs.GridFunction(self.codomain.fes) help_mu2=ngs.GridFunction(self.codomain.fes) help_ET.Set(tuple(map_ET)) help_EA.Set(tuple(map_EA)) help_mu2.Set(tuple(map_mu2)) return help_EA, help_ET, help_mu2 def _eval(self,argument): self._read_in(argument, self.gfu_a) #scaling factor scale is approximately the ratio between E_A and E_T times nu_A. This has to be fixed in order to obtain #a PDE which is linear in the material constants #elastic constants E_A=self.gfu_a.components[0] E_T=self.gfu_a.components[1] D=(1-self.nu_T-2*self.scale*self.nu_A)*(1+self.nu_T) lam1=(1-self.nu_T**2)*E_A/D lam2=(1-self.nu_A*self.scale)*E_T/D lam4=(self.nu_A+self.nu_A*self.nu_T)*E_T/D lam6=(self.nu_T+self.nu_A*self.scale)*E_T/D self.list_mat[0]=lam1 if self.dim==2: mu2=0.5*E_A/(1+self.nu_T) self.list_mat[15]=lam2 if self.dim==3: mu1=0.5*E_T/(1+self.nu_T) mu2=self.gfu_a.components[2] for i in self.list_lam6: self.list_mat[i]=lam6 for i in self.list_mu1: self.list_mat[i]=mu1 for i in self.list_lam2: self.list_mat[i]=lam2 #definition of material tensor self.list_mat[0]=lam1 for i in self.list_lam4: self.list_mat[i]=lam4 for i in self.list_mu2: self.list_mat[i]=mu2 #convert to gridfunction for further use self.gfu_tens.Set(tuple(self.list_mat)) #flat numpy-array from gridfunction return self.gfu_tens.vec.FV().NumPy().copy() def _adjoint(self,argument): self._read_in(argument, self.gfu_map) adj_ET=ngs.InnerProduct(self.gfu_map,self.help_ET) if self.dim==2: adj_EA=ngs.InnerProduct(self.gfu_map,self.help_EA+self.help_mu2) help2=(adj_EA,adj_ET) if self.dim==3: adj_EA=ngs.InnerProduct(self.gfu_map,self.help_EA) adj_mu2=ngs.InnerProduct(self.gfu_map, self.help_mu2) help2=(adj_EA,adj_ET,adj_mu2) self.gfu_a.Set(help2) return self.gfu_a.vec.FV().NumPy().copy()maps physical input parameters to subspace in material tensor space transverse isotropic elastice with isotropic plane 23
Ancestors
Inherited members
class ProjectToBoundary (domain, codomain=None)-
Expand source code
class ProjectToBoundary(NGSolveOperator): def __init__(self, domain, codomain=None): codomain = codomain or domain super().__init__(domain, codomain) self.linear=True self.bdr = codomain.bdr self.gfu_codomain = ngs.GridFunction(self.codomain.fes) self.gfu_domain = ngs.GridFunction(self.domain.fes) try: self.nr_bc = len(self.codomain.summands) except: self.nr_bc = 1 def _eval(self, x): if self.nr_bc == 1: array = [x] else: array = self.domain.split(x) toret = [] for i in range(self.nr_bc): self.gfu_domain.vec.FV().NumPy()[:] = array[i] self.gfu_codomain.Set(self.gfu_domain, definedon=self.codomain.fes.mesh.Boundaries(self.bdr)) toret.append(self.gfu_codomain.vec.FV().NumPy().copy()) return np.array(toret).flatten() def _adjoint(self, g): toret = [] if self.nr_bc == 1: g_tuple = [g] else: g_tuple = self.codomain.split(g) for i in range(self.nr_bc): self.gfu_codomain.vec.FV().NumPy()[:] = g_tuple[i] self.gfu_domain.Set(self.gfu_codomain, definedon=self.codomain.fes.mesh.Boundaries(self.bdr)) toret.append(self.gfu_domain.vec.FV().NumPy().copy()) return np.array(toret).flatten()Ancestors
Inherited members
class EIT (domain, g, codomain=None, alpha=0.01)-
Expand source code
class EIT(NGSolveOperator): r"""Electrical Impedance Tomography Problem PDE: \[ -\textrm{div}(s \nabla u)+\alpha u=0 \;\text{ in } \Omega \] \[ s \frac{\textrm{d}u}{\textrm{d}n} = g \;\text{ on } \partial\Omega \] Evaluate: \(F\colon s \mapsto \mathrm{tr}(u)\) Derivative: \[ -\textrm{div}(s \nabla v)+\alpha v=\textrm{div}(h \nabla u) (=:f) \] \[ s \frac{\textrm{d}v}{\textrm{d}n} = 0 +(-h\frac{\textrm{d}u}{\textrm{d}n} \;\text{ [second term often omitted] } \] Der: \(F'[s]\colon h \mapsto \textrm{tr}(v)\) Adjoint: \[ -\textrm{div}(s \nabla w)+\alpha w=0 \] \[ s \frac{\textrm{d}w}{\textrm{d}n} = q \] Adj: \(F'[s]^*\colon q \mapsto -\nabla(u) \nabla(w)\) Proof: \[(F'h, q)=\int_{\partial\Omega} [\textrm{tr}(v) q] \] \[= \int_{\partial\Omega} [\textrm{tr}(v) s \frac{\textrm{d}w}{\textrm{d}n}] \] \[= \int_{\Omega} [\textrm{div}(v s \nabla w )]\] Note \(\textrm{div}(s \nabla w) = \alpha*w\), thus above equation shows: \[(F'h, q) = (s \nabla v, \nabla w)+\alpha (v, w) \] \[= \int_\Omega [\textrm{div}( s \nabla v w)] +(-\textrm{div} (s \nabla v)), w)+\alpha (v, w)\] \[= \int_{\partial\Omega} [s dv/dn \textrm{tr}(w)]+(f, w)\] \[= (f, w)-\int_{\partial\Omega} [\textrm{tr}(w) h \frac{\textrm{d}u}{\textrm{d}n}]\] \[= (h, -\nabla u \nabla w) + \int_\Omega [\textrm{div}(h \nabla u w)]-\int_{\partial\Omega} [\textrm{tr}(w) h \frac{\textrm{d}u}{\textrm{d}n}] \] The last two terms are the same! It follows: \((F'h, q) = (h, -\nabla u \nabla w)\). Hence: Adjoint: \(q \mapsto -\nabla u \nabla w\) """ def __init__(self, domain, g, codomain=None, alpha=0.01): codomain = codomain or domain #Need to know the boundary to calculate Neumann bdr condition assert codomain.bdr is not None super().__init__(domain, codomain) self.g = g self.nr_bc = len(self.g) self.fes_domain = domain.fes self.fes_codomain = codomain.fes #FES and Grid Function for reading in values self.fes_in = ngs.H1(self.fes_codomain.mesh, order=1) self.gfu_in = ngs.GridFunction(self.fes_in) # grid functions for later use self.gfu_eval = ngs.GridFunction(self.fes_codomain) # solution, return value of _eval self.gfu_deriv = ngs.GridFunction(self.fes_codomain) # grid function return value of derivative self.gfu_adjoint = ngs.GridFunction(self.fes_domain) #grid function return value of adjoint self.gfu_bf = ngs.GridFunction(self.fes_codomain) # grid function for defining integrator (bilinearform) self.gfu_lf = ngs.GridFunction(self.fes_codomain) # grid function for defining right hand side (linearform), f self.gfu_b = ngs.GridFunction(self.fes_codomain) self.gfu_inner_adjoint = ngs.GridFunction(self.fes_codomain) # grid function for inner computations in adjoint self.Number = ngs.NumberSpace(self.fes_codomain.mesh) #r, s = self.Number.TnT() u, v = self.fes_codomain.TnT() # Define Bilinearform, will be assembled later self.a = ngs.BilinearForm(self.fes_codomain, symmetric=True) self.a += (ngs.grad(u) * ngs.grad(v) * self.gfu_bf+alpha*u*v) * ngs.dx #Additional condition: The integral along the boundary vanishes #self.a += ngs.SymbolicBFI(u * s + v * r, definedon=self.fes_codomain.mesh.Boundaries("cyc")) #self.fes1 = ngs.H1(self.fes_codomain.mesh, order=4, definedon=self.fes_codomain.mesh.Boundaries("cyc")) # Define Linearform for evaluation, will be assembled later self.b = ngs.LinearForm(self.fes_codomain) self.b += self.gfu_b*v*ngs.ds(codomain.bdr) # Define Linearform for derivative, will be assembled later self.f_deriv = ngs.LinearForm(self.fes_codomain) self.f_deriv += -self.gfu_lf * ngs.grad(self.gfu_eval) * ngs.grad(v) * ngs.dx # Initialize preconditioner for solving the Dirichlet problems by ngs.solvers.BVP self.prec = ngs.Preconditioner(self.a, 'direct') #Weak formulation: #0=int_Omega [-div(s grad u) v + alpha u v]=-int_dOmega [s du/dn trace(v)]+int_Omega [s grad u grad v + alpha u v] #Hence: int_Omega [s grad u grad v + alpha u v] = int_dOmega [g trace(v)] #Left term: Bilinearform self.a #Righ term: Linearform self.b def _eval(self, diff, differentiate=False, adjoint_derivative=False): # Assemble Bilinearform self._read_in(diff, self.gfu_bf) self.a.Assemble() # Assemble Linearform, boundary term toret = [] for i in range(self.nr_bc): self.gfu_b.Set(self.g[i]) self.b.Assemble() # Solve system if i == 0: self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_eval, prec=self.prec) else: self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_eval, prec=self.prec) toret.append(self.gfu_eval.vec.FV().NumPy()[:].copy()) return np.array(toret).flatten() #Weak Formulation: #0 = int_Omega [-div(s grad v) w + alpha v w]-int_Omega [div (h grad u) w] #=-int_dOmega [s dv/dn trace(w)] + int_Omega [s grad v grad w + alpha v w]-int_dOmega [h du/dn trace(w)]+int_Omega [h grad u grad w] #=int_Omega [s grad v grad w + alpha v w]+int_Omega [h grad u grad w] #Hence: int_Omega [s grad v grad w + alpha v w] = int_Omega [-h grad u grad w] #Left Term: Bilinearform self.a, already defined in _eval #Right Term: Linearform f_deriv def _derivative(self, h, **kwargs): # Bilinearform already defined from _eval # Assemble Linearform toret = [] for i in range(self.nr_bc): self._read_in(h[i], self.gfu_lf) self.f_deriv.Assemble() self.gfu_deriv.Set(0) self._solve_dirichlet_problem(bf=self.a, lf=self.f_deriv, gf=self.gfu_deriv, prec=self.prec) toret.append(self.gfu_deriv.vec.FV().NumPy()[:].copy()) return np.array(toret).flatten() #Same problem as in _eval def _adjoint(self, argument): # Bilinearform already defined from _eval # Definition of Linearform # But it only needs to be defined on boundary if self.nr_bc==1: argument_tuple = [argument] else: argument_tuple = self.codomain.split(argument) toret = np.zeros(np.size(self.gfu_adjoint.vec.FV().NumPy())) for i in range(self.nr_bc): self.gfu_b.vec.FV().NumPy()[:] = argument_tuple[i] self.b.Assemble() self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_inner_adjoint, prec=self.prec) self.gfu_adjoint.Set(-ngs.grad(self.gfu_inner_adjoint) * ngs.grad(self.gfu_eval)) toret += self.gfu_adjoint.vec.FV().NumPy().copy() return toretElectrical Impedance Tomography Problem
PDE: -\textrm{div}(s \nabla u)+\alpha u=0 \;\text{ in } \Omega s \frac{\textrm{d}u}{\textrm{d}n} = g \;\text{ on } \partial\Omega Evaluate: F\colon s \mapsto \mathrm{tr}(u) Derivative: -\textrm{div}(s \nabla v)+\alpha v=\textrm{div}(h \nabla u) (=:f) s \frac{\textrm{d}v}{\textrm{d}n} = 0 +(-h\frac{\textrm{d}u}{\textrm{d}n} \;\text{ [second term often omitted] }
Der: F'[s]\colon h \mapsto \textrm{tr}(v)
Adjoint: -\textrm{div}(s \nabla w)+\alpha w=0 s \frac{\textrm{d}w}{\textrm{d}n} = q
Adj: F'[s]^*\colon q \mapsto -\nabla(u) \nabla(w)
Proof: (F'h, q)=\int_{\partial\Omega} [\textrm{tr}(v) q] = \int_{\partial\Omega} [\textrm{tr}(v) s \frac{\textrm{d}w}{\textrm{d}n}] = \int_{\Omega} [\textrm{div}(v s \nabla w )] Note \textrm{div}(s \nabla w) = \alpha*w, thus above equation shows: (F'h, q) = (s \nabla v, \nabla w)+\alpha (v, w) = \int_\Omega [\textrm{div}( s \nabla v w)] +(-\textrm{div} (s \nabla v)), w)+\alpha (v, w) = \int_{\partial\Omega} [s dv/dn \textrm{tr}(w)]+(f, w) = (f, w)-\int_{\partial\Omega} [\textrm{tr}(w) h \frac{\textrm{d}u}{\textrm{d}n}] = (h, -\nabla u \nabla w) + \int_\Omega [\textrm{div}(h \nabla u w)]-\int_{\partial\Omega} [\textrm{tr}(w) h \frac{\textrm{d}u}{\textrm{d}n}] The last two terms are the same! It follows: (F'h, q) = (h, -\nabla u \nabla w). Hence: Adjoint: q \mapsto -\nabla u \nabla w
Ancestors
Inherited members
class ReactionNeumann (domain, g, codomain=None)-
Expand source code
class ReactionNeumann(NGSolveOperator): r""" Estimation of the reaction coefficient from boundary value measurements PDE: -div(grad(u)) + s*u = 0 in Omega du/dn = g on dOmega Evaluate: F: s \mapsto trace(u) Derivative: -div(grad(v))+s*v = -h*u (=:f) dv/dn = 0 Der: F'[s]: h \mapsto trace(v) Adjoint: -div(grad(w))+s*w = 0 dw/dn = q Adj: F'[s]^*: q \mapsto -u*w proof: (F'h, q) = int_dOmega [trace(v) q] = int_dOmega [trace(v) dw/dn] = int_Omega [div(v grad w)] = int_Omega [grad v grad w] + int_Omega [v div( grad w)] = int_Omega [div(w grad v)] - int_Omega [div(grad v) w] + int_Omega [v div (grad w)] = int_dOmega [trace(w) dv/dn] - int_Omega [h u w] - int_Omega [s v w] + int_Omega [v s w] Note that dv/dn=0 on dOmega. Hence: (F'h, q) = -int_Omega[h u w] = (h, -u w) """ def __init__(self, domain, g, codomain=None): codomain = codomain or domain #Need to know the boundary to calculate Neumann bdr condition assert codomain.bdr is not None super().__init__(domain, codomain) self.g = g self.nr_bc = len(self.g) self.fes_domain = domain.fes self.fes_codomain = codomain.fes # grid functions for later use self.gfu_eval = ngs.GridFunction(self.fes_codomain) # solution, return value of _eval self.gfu_deriv = ngs.GridFunction(self.fes_codomain) # grid function: return value of derivative self.gfu_adjoint = ngs.GridFunction(self.fes_domain) # grid function: return value of adjoint self.gfu_bf = ngs.GridFunction(self.fes_codomain) # grid function for defining integrator of bilinearform self.gfu_lf = ngs.GridFunction(self.fes_domain) # grid function for defining linearform self.gfu_b = ngs.GridFunction(self.fes_codomain) # grid function for defining the boundary term self.gfu_inner_adjoint = ngs.GridFunction(self.fes_codomain) # grid function for inner computation in adjoint #Test and Trial Function u, v = self.fes_codomain.TnT() # Define Bilinearform, will be assembled later self.a = ngs.BilinearForm(self.fes_codomain, symmetric=True) self.a += (ngs.grad(u) * ngs.grad(v) + u * v * self.gfu_bf) * ngs.dx # Boundary term self.b = ngs.LinearForm(self.fes_codomain) self.b += -self.gfu_b * v.Trace() * ngs.ds(codomain.bdr) # Linearform (only appears in derivative) self.f_deriv = ngs.LinearForm(self.fes_codomain) self.f_deriv += -self.gfu_lf * self.gfu_eval * v * ngs.dx # Initialize preconditioner for solving the Dirichlet problems by ngs.solvers.BVP self.prec = ngs.Preconditioner(self.a, 'direct') def _eval(self, diff, differentiate=False, adjoint_derivative=False): # Assemble Bilinearform self._read_in(diff, self.gfu_bf) self.a.Assemble() # Assemble Linearform of boundary term toret = [] for i in range(self.nr_bc): self.gfu_b.Set(self.g[i]) self.b.Assemble() # Solve system if i == 0: self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_eval, prec=self.prec) else: self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_eval, prec=self.prec) toret.append(self.gfu_eval.vec.FV().NumPy()[:].copy()) return np.array(toret).flatten() def _derivative(self, h): # Bilinearform already defined from _eval # Assemble Linearform of derivative toret = [] for i in range(self.nr_bc): self._read_in(h, self.gfu_lf) self.f_deriv.Assemble() # Solve system self._solve_dirichlet_problem(bf=self.a, lf=self.f_deriv, gf=self.gfu_deriv, prec=self.prec) toret.append(self.gfu_deriv.vec.FV().NumPy()[:].copy()) return np.array(toret).flatten() def _adjoint(self, argument): # Bilinearform already defined from _eval # Definition of Linearform # But it only needs to be defined on boundary if self.nr_bc==1: argument_tuple = [argument] else: argument_tuple = self.codomain.split(argument) toret = np.zeros(np.size(self.gfu_adjoint.vec.FV().NumPy())) for i in range(self.nr_bc): self.gfu_b.vec.FV().NumPy()[:] = argument_tuple[i] self.b.Assemble() # Solve system self._solve_dirichlet_problem(bf=self.a, lf=self.b, gf=self.gfu_inner_adjoint, prec=self.prec) self.gfu_adjoint.Set(self.gfu_inner_adjoint * self.gfu_eval) toret+=self.gfu_adjoint.vec.FV().NumPy().copy() return toretEstimation of the reaction coefficient from boundary value measurements
PDE: -div(grad(u)) + s*u = 0 in Omega du/dn = g on dOmega
Evaluate: F: s \mapsto trace(u)
Derivative
-div(grad(v))+sv = -hu (=:f) dv/dn = 0
Der: F'[s]: h \mapsto trace(v)
Adjoint: -div(grad(w))+sw = 0 dw/dn = q Adj: F'[s]^: q \mapsto -u*w
proof: (F'h, q) = int_dOmega [trace(v) q] = int_dOmega [trace(v) dw/dn] = int_Omega [div(v grad w)] = int_Omega [grad v grad w] + int_Omega [v div( grad w)] = int_Omega [div(w grad v)] - int_Omega [div(grad v) w] + int_Omega [v div (grad w)] = int_dOmega [trace(w) dv/dn] - int_Omega [h u w] - int_Omega [s v w] + int_Omega [v s w] Note that dv/dn=0 on dOmega. Hence: (F'h, q) = -int_Omega[h u w] = (h, -u w)
Ancestors
Inherited members