Hot answers tagged

12

You can interpolate the solution onto a finer mesh and then plot it: from dolfin import * coarse_mesh = UnitSquareMesh(2, 2) fine_mesh = refine(refine(refine(coarse_mesh))) P2_coarse = FunctionSpace(coarse_mesh, "CG", 2) P1_fine = FunctionSpace(fine_mesh, "CG", 1) f = interpolate(Expression("sin(pi*x[0])*sin(pi*x[1])"), P2_coarse) g = interpolate(f, ...


11

I suggest you start by looking at the FEniCS Navier-Stokes demo which is documented here: http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/navier-stokes/python/documentation.html For your specific test case, you might want to look at the NSBENCH set of test problem (which are described in Chapter 21 of the FEniCS Book). The code for that ...


10

To implement your problem in FEniCS, you have to replace the integrals in terms of boundaries by integrals in terms of edges. This introduces jumps/averages in the test functions, which you entirely miss in your implementation. Hence, the system is not invertible and your solution does not look right. Equation (3.3) in Arnold et. al. 2002 gives you a tool to ...


10

I've been developing a lightweight finite element library in Python 2.7 harnessing the power of NumPy arrays and SciPy sparse matrices. The general idea is that given a mesh and a finite element, you have more-or-less one-to-one correspondence between the bilinear form and a (sparse) matrix. The user can then use the resulting matrix as he or she sees fit. ...


8

First, a general point: you cannot prescribe arbitrary boundary conditions for a partial differential operator and expect that the partial differential equation (which always includes both operator and boundary conditions) is well-posed, i.e., admits a unique solution that depends continuously on the data -- all of which is a necessary condition for actually ...


8

That depends how did you generated the file. For example your mesh generator can produce some cell or facet markers according to its input. Then dolfin-convert script may or may not succeed converting all these markers to XML. They're then stored within XML mesh file or in separate XML. You can then use them as FacetFunction # if stored within mesh ...


8

x = V.cell().x then use x[0] and x[1] as $x$ and $y$ respectively. Right-hand side of LinearVariationalProblem needs to be rank 1 form - this is expressed by dependendnce on TestFunction and independence on TrialFunction.


8

I've been working a bit on adaptive refinement to do the job (see the code below). The scaling of the error indicator with total mesh size and total variation of the mesh function is not perfect, but you could fit that to your needs. The images below are for testcase #4. The number of cells increases from 200 to about 24,000, which may be a bit over the top, ...


8

As one of the library's authors, I would of course love for deal.II to come out on top with this comparison. But I suspect it may not, and the answer lies in a factor you omit from your comparison: how long it actually takes to implement your code. Few people in academia with the skills to implement a FEM code from scratch spend more time solving PDEs than ...


7

Use dolfin-convert from a terminal window.


7

Part of your problem is your choice of boundary conditions. This paper (which I wrote) develops the criteria for imposing a pressure boundary condition on the incompressible Navier-Stokes equations. The short version is that you need to impose the normal component of the normal traction on the inlet and outlet to be the pressures you want on each surface, ...


7

You can use a Nitsche-type method for this. See the following reference: J. Freund, R. Stenberg. On weakly imposed boundary conditions for second order problems. Proceedings of the Ninth Int. Conf. Finite Elements in Fluids, Venice 1995. M. Morandi Cecchi et al., Eds. pp. 327-336. I have implemented this a while ago in some simple FEniCS-Code to deal with ...


7

Doing u, p = U.split() gives you a view into U, i.e. u, p and U share degree of freedom data. Doing u, p = U.split(deepcopy=True) gives a u and a p that have their own degree-of-freedom vectors. Use of shallow and deep sub-functions is demonstrated in the DOLFIN demo demo/undocumented/stokes-taylor-hood (from both C++ and Python). The syntax is a bit ...


6

Your sample code will not work in parallel because you are creating the same vertices and cells on all processes. If you add mesh entities only on process 0 and then call MeshPartitioning.build_distributed_mesh(mesh) on all processes, your sample code runs fine.


6

I've used this exact workflow for solving the Navier-Stokes equations in FEniCS. Almost certainly the hardest part in the process will be implementing the numerical solution method; there is no out-of-the-box support for that in FEniCS (although of course all differential operators you need are available). For reference, check out nsbench, a Navier-Stokes ...


5

You can always use an iterator to access the entities connected to a cell, like the edges in your case. Just as you suggest, you first need to create a dolfin::Cell from the ufc::cell. Then do something like this: void eval(Array<double>& v, const Array<double>& x, ...


5

Looks like you're using the Python interface, in which case u is both a DOLFIN Function object and a UFL object (it inherits from both). The type information you're seeing is related to UFL. You can evaluate the function u at an arbitrary point: f = u((0.5, 0.5, 0.5)) if DOLFIN is configured with CGAL. I don't know how you could check that the solution ...


5

That's actually very easy to do in Dolfin: from dolfin import * # define mesh, function space (piecewise linear) mesh = UnitSquareMesh(64,64) V = FunctionSpace(mesh,'CG',1) # inhomogeneous boundary conditions (otherwise the solution is trivial) bc = DirichletBC(V, Constant(1.0), lambda x,on_boundary: on_boundary) # define bi(non)linear form # note that ...


5

This is I think rotating cylinder and Karman vortex street behind it. You can also uncomment SUPG/PSPG stabilization. But I suspect that mesh size h computed as 2*circumradius is not good quantity and prefer cell diameter computed this way https://scicomp.stackexchange.com/a/7181/4254 from dolfin import * import time #set_log_level(PROGRESS) parameters['...


5

As the error message suggest the DirichletBC does not hit any mesh entities with corresponding dofs. You need to examine your subdomains. One way to debug these are to apply the DirichletBC to a Function and plot the result: def left_boundary(x, on_boundary): r = math.sqrt(x[0] * x[0] + x[1] * x[1]) return r < 0.5 leftPlate = DirichletBC(V, u_L, ...


5

Since the problem is treated in subdomains-poisson for the case, where the mesh is aligned with the interface, I assume in the following that in your case the jump can occur anywhere. I think the first step here is not to think about implementation FEniCS, but rather to find an adequate variational form for this problem, where the interface conditions are ...


5

Take a look at the CSG (Constructive Solid Geometry) demos in DOLFIN. You should be able to generate a rectangle with two holes with something like mesh = Mesh(Rectangle(...) - Ellipse(...) - Ellipse(...), resolution)


5

Make sure that the mesh is aligned with x = 0. Then mark the facets along that line and compute the flux using an "interior facet integral". You can do this by supplying a facet function markers that (e.g.) marks those facets by 1 and the remaining facets by 0. Then compute the flux by line_segment = AutoSubDomain(lambda x: near(x[0], 0)) markers = ...


5

You need to evaluate the Function at the given point. Assuming phi1 is a scalar Function, try: class Bulk(SubDomain): def inside(self, x, on_boundary): return phi1(x) >= 0.1 and similar for the other SubDomains


4

The following compiles: from dolfin import * from random import random lmbda = Constant(1.0) mu = Constant(1.0) def sigma(v): return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d) mesh = UnitCubeMesh(8,8,8) V = VectorFunctionSpace(mesh, "CG", 1) # Create random displacement vector displacement = Function(V) for i in xrange(V.dim()...


4

You can define the Mini element in this alternative manner if you want a mixed enriched element rather than an enriched mixed element: U = FunctionSpace(mesh, "Lagrange", 1) B = FunctionSpace(mesh, "Bubble", 3) M = U + B Mini_h = MixedFunctionSpace([M, M])


4

There will be non-negligible overhead when moving between meshes (some current developments will make this faster than it is now). You could first interpolate, then to evaluate vertex values efficiently use the function Function::compute_vertex_values. For spaces that are not continuous it can yield some unexpected results.


4

I don't know about these things in FEniCS. For openFOAM there are libraries for multiphase flow. The implementation of population balance equations to track particle size distributions in openFOAM has been done and discussed in this paper. The question 'openFOAM or FEniCS' can maybe decided by deciding on whether you want to use finite volumes or finite ...


4

You can't expect that solution to your altered problem would be a solution to Poisson problem because you need to change the problem somehow to make it well-posed. One could guess that possible problem formulation is to minimize $$F(u, \lambda) = \int_\Omega \frac{1}{2}|\nabla u|^2 \;\mathrm{d}x - \int_\Omega f u \;\mathrm{d}x - \int_{\Gamma_\mathrm{N}} g ...


4

Vatiational formulation of Poisson problem with $\varepsilon \in L^\infty$ Both interfacial conditions are incorporated in following variational problem. Given $\Omega$ Lipschitz domain $\partial\Omega = \Gamma_\mathrm{D} \cup \Gamma_\mathrm{N}$ $V=\{v\in H^1(\Omega); v|_{\Gamma_\mathrm{D}}=0 \} $ $\varepsilon\in L^\infty(\Omega)$ $g \in L^2(\Gamma_\...


Only top voted, non community-wiki answers of a minimum length are eligible