Extracting results from a boundary


175
views
0
5 months ago by
Elch  
Hello Community,

I'm new to the forum and i'm quite sure that my problem was allready solved, yet I'm to stupid to find the solution.

So let's start, I have a mesh generated by gmsh with several subdomains/boundaries. After importing the mesh, i'll apply Dirilichet BC (Temperatures for this case) on all boundaries and solve the heat equation.
(Problem is static heat conduction)
After calculating the Temperaturedistribution, I'll further calculate the heat flux over the domain.

Now what really interests me is the heat flux on my boundaries, but I find no way to extract these from the results...

class bc_north(Expression):
    def set_temp_north(self, par):
        self.par = par
    def eval(self, value, x):
        value[0]= x[0]*x[0]*self.par[0] + x[0]*self.par[1]+self.par[2]

class bc_south(Expression):
    def set_temp_south(self, par):
        self.par = par
    def eval(self, value, x):
        value[0]=x[0]*x[0]*self.par[0] + x[0]*self.par[1]+self.par[2]

class bc_east(Expression):
    def set_temp_east(self, par):
        self.par = par
    def eval(self, value, x):
        value[0]= x[1]*x[1]*x[1]*self.par[0]+x[1]*x[1]*self.par[1]+x[1]*self.par[2]+self.par[3]


def sim_A(north_par,south_par,east_par,west_par):
    ###loading mesh
    mesh = Mesh('DiskA.xml')
    boundaries = MeshFunction("size_t", mesh,'DiskA_facet_region.xml')
    subdomains = MeshFunction("size_t", mesh,'DiskA_physical_region.xml')
    element = FiniteElement("CG", tetrahedron,2)
    F = FunctionSpace(mesh, element) 
    V = VectorFunctionSpace(mesh, element, degree=1,dim=3)
    u = TrialFunction(F)
    v = TestFunction(F)
    r = Expression('x[1]', degree=1)
    #settin bc´s
    southtemp = bc_south(degree=2)
    southtemp.set_temp_south(south_par)
    bc1 = DirichletBC(F, southtemp, boundaries, 1)
    northtemp = bc_north(degree=2)
    northtemp.set_temp_north(north_par)
    bc2 = DirichletBC(F, northtemp, boundaries, 2)
    easttemp = bc_east(degree=2)
    easttemp.set_temp_east(east_par)
    bc3 = DirichletBC(F, easttemp, boundaries, 4)
    westtemp = bc_east(degree=2)
    westtemp.set_temp_east(west_par)
    bc4 =DirichletBC(F, westtemp, boundaries, 5)
    conductivity = Constant(0.015)
    #defining the problem and solving
    a = conductivity*(Dx(u,0)*Dx(v,0)+Dx(u,1)*Dx(v,1))*r*dx()
    f=Constant(0)
    L = f*v*dx
    u = Function(F)
    solve(a==L, u, [bc1, bc2, bc3, bc4])
    #calculating the heat flux
    grad_u = project(conductivity*-grad(u), V)
#    this one does the trick, but I dont want the flux integrated
#    n = FacetNormal(mesh)
#    ds = Measure("ds")[boundaries]
#    flux = assemble(dot(grad_u,n)*ds(4))
    vtkfile = File('A_Temp.pvd')
    vtkfile << u
    vtkfile = File('A_flux.pvd')
    vtkfile << grad_u
​


I hope my code is not to long for the forum and i'll take any advice for a cleaner code thankfully.
The part with ds(4) is quite the right direction but integrates the flux.
I need more like a vector, my goal is to plot the flux on the boundaries over x[1] (or x[0] for north and south).

I found some crazy stuff with submeshs and mapping dofmaps from mesh to submesh and vice versa in the old forums, but wasnt able to reproduce. I think at this point some syntax changes hit me.

Just for the protocoll, I'm not a native speaker, but i hope the english is not all that bad.

I would be thankfull for every advice i might get.

Greetings

Elch

Community: FEniCS Project
Have you found a solution in the meantime?
written 5 weeks ago by polyeder  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »