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

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

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

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

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

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

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

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

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

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

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

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_\...


4

So the problem is that FEniCS doesn't recognize that in 1D, u and grad(u) have the same shape. (You could file an issue about this on the Dolfin issue tracker.) If you replace grad by div, your code works. (The general form of the inviscid Burgers equation is $\mathop{\mathrm{div}}(u^2) = f$.)


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'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

Use the highest order you need for all terms. The major expense is evaluating the values and gradients of the shape functions at the quadrature points, which you do once and then reuse as often as is necessary. If you use multiple different quadrature formulas, you end up with more quadrature points at which you need to evaluate shape functions, for a larger ...


4

I would try pushing this form V = VectorFucntioSpace(mesh, 'CG', 1, dim=2) u = Function(V) v = TestFunction(V) sigma = Constant(0.1) F = inner(grad(u), grad(v))*dx \ + (inner(u, u)-1.0)*((u[0]+sigma*u[1])*v[0] + (u[1]-sigma*u[0])*v[1])*dx into nonlinear-poisson demo. But note that this would yield trivial solution if you start from zero field with zero ...


4

You might have to reassemble if your problem is non-linear and your method at a future step incorporates the solution in the formation of the matrix. If you are doing Picard iteration rather than Netwon-Raphson, then you should only have to reform the right-hand-side vector. I don't know enough about FEciCS and COMSOL to say what they do, but I suspect, for ...


4

You seem to be confused about which equation to solve. You have two: (i) the flow equations, (ii) the equations for your particle property $D$. The finite element method is suitable for part 1 of this. Or, if you really just have a pipe with laminar flow, then you actually know what the flow field is (namely, Poisseuille flow) and you don't need to solve ...


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