Extracting results from a boundary
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= x*x*self.par + x*self.par+self.par class bc_south(Expression): def set_temp_south(self, par): self.par = par def eval(self, value, x): value=x*x*self.par + x*self.par+self.par class bc_east(Expression): def set_temp_east(self, par): self.par = par def eval(self, value, x): value= x*x*x*self.par+x*x*self.par+x*self.par+self.par 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', 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 (or x 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.