Hot answers tagged

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

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


4

I think you have some confusion. PETSc is not in the same league as Fenics, Libmesh, Moose etc. In fact, all of these (heavyweight) packages use PETSc for linear algebra. IMHO PETSc is as lightweight as you can get. It just requires C/Fortran compilers and Python (used only for configuration) and you can build the library in under 5 minutes on your laptop. ...


4

Specific answers to this question are probably time-limited. However, the following general approach (from the great Eric S. Raymond) works very well: Rule of Modularity: Write simple parts connected by clean interfaces. Rule of Clarity: Clarity is better than cleverness. Rule of Composition: Design programs to be connected to other programs. Rule of ...


4

Unfortunately there's no tool for this. You can run each on a variety of input sizes to establish the computational complexity they appear to have, i.e. the $f$ in the $O(f(n))$ that characterizes each code. This can point you towards what underlying algorithms each is using and verify or not that something that should be $O(n)$ is actually implemented to ...


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


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

The SEM is a FEM! It's almost like all these different names are designed to confuse the newcomer. I will speak primarily about the most popular form which uses a tensor product Lagrange basis with nodes at the Legendre-Gauss-Lobatto (LGL) points. Once you have chosen that basis and integrate with the LGL quadrature using the same nodes. You have a method ...


4

The main advantage is that it reduces the Runge phenomenon and leads to faster convergence rates. It also presents less numerical dispersion and need less nodes per wavelength (see 1 and 2). So, I would say that you would prefer the method for wave propagation scenarios. Regarding software that includes SEM, I am aware of the following: FSELib: Matlab ...


4

I can recommend nutils. nutils meets at least a few your "light-weight" requirements. it is pure python and easy to install since it only depends on standard Python libraries numpy, scipy, and matplotlib and, thus, it is well suited for interoperations. At least the developers claim that "Exposed objects are of native python type or allow for easy ...


3

My suggestion is to just write your own code using PETSc instead of using an existing FE library. Parallel assembly/solve is the most complicated part of an FE code and PETSc takes care of it. The rest of it is simple anyway, maybe a few hundred lines of straightforward C/Fortran. Besides PETSc now has DMPLEX for managing meshes (though it is still in ...


3

To find the weaks forms we use the identity $$ \nabla\cdot(\rho\,\nabla\phi) = \rho\,\nabla^2\phi + \nabla\rho \cdot \nabla\phi $$ and its vector counterpart $$ \nabla\cdot(\varphi\,\mathbf{v}) = \varphi\,(\nabla\cdot\mathbf{v}) + (\nabla\varphi) \cdot \mathbf{v} \,. $$ Then the first equation can be written as $$ \partial_{t}\rho+2\,\nabla\cdot(\rho\,...


3

FEM is for solving boundary value problems. What you have here is an initial value problem (assuming $\alpha$ and $U_z$ are constant or functions of $z$), the same as solving a time-dependent ODE, except here $z$ is your time-like coordinate. The appropriate way to solve this kind of problem is with time-integration methods, e.g. Runge-Kutta schemes. In FEM ...


2

Refer to dolfin-convert script installed within FEniCS, and its manual page, $ man dolfin-convert for supported mesh formats. Then pick some software generating those formats. Note that cell functions (useful for specifying permitivities in your case) and/or facet functions (useful for boundary and interfacial conditions) are not fully supported by dolfin-...


2

I have also had this problem and spent a lot of time on various forums and I have finally come up with a good solution. This solution will allow you to specify a boundary condition on any whole subdomain within your geometry, and each subdomain can have its own separate boundary condition. I am posting this solution here as the old official Fenics QA forum ...


2

meshio (a small project of mine) has as command-line tool that can do it: meshio-convert cylinder6.msh out.xml


2

There is currently no documented 'FEniCS way' to do this. However, since FEniCS is a pretty standard finite element code behind all the UFL and code generation magic, you can implement things like transfer operators by yourself. The only difficulty is that they have no built-in mechanism to deal with the inter-level mappings in hierarchically refined meshes ...


2

If you look at the documentation,1 $u$ is described as solving the Poisson equation $(\nabla u,\nabla v) = (f,v)$ for $v\in H^1(\Omega)$. This means that it is unlikely to be a velocity field (it's scalar, after all). I would interpret the documentation insofar as only the geometric domain is meant to model a blood vessel, while the PDE is only a very ...


2

I agree with the suggestion of starting with a simple problem and with the elastic solution. Probably the simplest wave problem is the 1D, infinite bar/string. The analytical solution to this problem is well-known, e.g. http://mathworld.wolfram.com/dAlembertsSolution.html You can model this with a strip of elements in the x-direction. You want zero-stress ...


2

I'd suggest to try it on your own. Do an expansion of your wavefunction in terms of spherical harmonics, $$ \psi(\mathbf r) \ = \ \sum_{\ell} R_\ell(r,t) \, Y_{\ell 0} (\theta,\phi)\,. $$ Note that I've set the index $m$ in $Y_{\ell m}$ to zero, in order to account for the symmetry of your Hamiltonian with respect to rotations around the $x-y$ plane. This ...


2

FEniCS users have solved this problem before, but keep in mind that FEniCS does not natively support complex numbers right now in its code. Therefore you have to make a workaround. See: https://fenicsproject.org/qa/9209/how-to-use-complex-numbers-iterative-solvers https://fenicsproject.org/qa/10671/mixed-function-spaces-for-complex-valued-problems-in-2016-0


2

TL;DR: Let $k$ be the diffusion coefficient, $\theta$ the minimum angle between any two edges of the mesh, $d$ the space dimension, and $p$ the polynomial degree of the finite element basis you're using. Then taking $$\alpha \ge \frac{2 \cdot d \cdot (p + d - 1)}{\cos\theta\cdot \tan\theta/2} \cdot \frac{\max k}{\min k}$$ should guarantee that the ...


2

Since GMSH was mentioned in the comments as a possible meshing tool, the common way to deal with it there is to create a Physical Line. The simplest .geo (modeling and meshing a square): cl=0.5; Point(1) = {1.,1.,0.,cl}; Point(2) = {-1.,1.,0.,cl}; Point(3) = {-1.,-1.,0.,cl}; Point(4) = {1.,-1.,0.,cl}; Line(1) = {1,2}; Line(2) = {2,3}; Line(3) = {3,4}; ...


2

The different specialized VertexFunction, EdgeFunction, FaceFunction, FacetFunction, CellFunction have indeed been deprecated and (at least) as of version 2018.2 removed. Instead, they have been merged into a general purpose MeshFunction which takes the kind of mesh object they're defined on as a (mandatory) parameter called the "topological dimension". To ...


2

The following approach should work. A = assemble(a) bc.apply(A) solver = LUSolver(A) solver.parameters['reuse_factorization'] = True # maybe not needed Now you can solve for any right hand side f b = assemble(f*v*dx) bc.apply(b) u = Function(V) solver.solve(u.vector(),b) This avoids assembly and LU factorization in every solve. Only LU solve is required ...


1

If you're using the latest version of FEniCS, the meshing tools are deprecated as you've noticed. These tools are now under mshr module. You can install it via terminal: conda install -c conda-forge mshr Additionally, if you'd like to convert meshes, the dolfin-convert became deprecated as well. The functions are now bundled in meshio module. Similarly: ...


1

Ok I will answer my own question: The problem is in the line solve(L==a,A,bc) which needs to be replaced by solve(a==L,A,bc). The two versions seem to non interchangeable. Doing so will result in the magnetic field B=project(curl(A),V) which is shown in the picture below:


1

Regarding the question of how to check whether $|v|<g$ (which I would recommend to split into $v(x)<g(x)$ and $v(x)>-g(x)$, since $|v|$ will in general not be a piecewise polynomial), this is rather easy for piecewise constant functions and piecewise linear functions (for which the maximum and minimum are always attained in the nodes) since there it ...


1

You could try the FEniCS book or the FEniCS Tutorial.


1

This is not really a program as much as a big repository of codes. John Burkardt from Florida State University maintains a rich collection of scripts in C++, Matlab and Fortran for a large range of problems. You can check out the FEM 2D library or go up a level and look at other codes listed there. These have good commenting on them and you can look at the ...


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