User login


You are here

EFG Matlab Routines

These used to be hosted at Northwestern, but the files were taken down some time ago. The original 1d and 2d Matlab routines for the element-free Galerkin method are now located at

These routines are described in detail in the paper

J. Dolbow and T. Belytschko (1998), "An Introduction to Programming the Meshless Element Free Galerkin Method," Archives of Computational Methods in Engineering, vol. 5, no. 3, pp. 207--242.

At the moment, the above page is a poor replica of the original, but at least the routines are there. I will add more when I find time.


Your programs are very interesting.They are useful for me.

Do you have the paper "J. Dolbow and T. Belytschko (1998), "An Introduction to Programming the Meshless Element Free Galerkin Method," Archives of Computational Methods in Engineering, vol. 5, no. 3" ?

Would you mind if you send it to me ? 

My email 

Thank you.Have a nice day!

I'm not sure I have a preprint I can post here or even an electronic version of the paper to send you.  If you send me your regular mail address (by email please), I will be happy to send you a standard reprint. 

Otherwise the article can be purchased online at this location.  

bqtinh's picture

Dear Prof. Dolbow

With your codes, this is appeared a error of "enorm" when you used (48x24):

ndivl = 48;ndivw = 24;

 The result:

enorm =


What do you think? I think that it is wrong.

Best regards!



Thanks for pointing this out Tinh.  I'm fairly certain that there is a bug in the code that was likely introduced when I transferred it from some old files (and had to retype many parts).  

I will look into this.  

bqtinh's picture

Dear Prof. Dolbow

Do you see any bug within your codes? I am looking forward to hear from you soon, best wishes




If we change the line as "ndivl = 48;ndivw = 24;", that means we set more nodals, so we also should set more GAUSS POINTS, we change the line:

"ndivlq = 10;ndivwq = 4;" as "ndivlq = 20;ndivwq = 8;", then the matrix m will not be singual.

I'm a graduate student come from China, JiangSu province. My subject is Topology Optimization of Continuum Structure Based on Meshless Method. Hope communicate with everyone. thanks!

my msn and email are all :

bqtinh's picture

Dear Prof. Dolbow

When you do for case 47x23 follow:

ndivl = 47;ndivw = 23;


It is appeared a error while running:


??? Reference to a cleared variable nt.

Error in ==> efg2d at 76
ltht = length(nt);


What is it? It is strange!

Best regards!






Hi Tinh,

Do you think the meshless method can be applicable to Fracture Mechanics of Steel material (elastoplastic mat.) ?





THE LINE "if(x(1,j)==Lb)" CHANGE AS "if(x(1,j)==((Lb/ndivl)*ndivl))"

THEN ok!

Dear Prof. Dolbow,

I'm Adam Othman, a 3rd year Mechanical Engineering Student in University of Bristol. Currently, I'm doing research project on EFG and need have to produce codes for a  linear elastic 2D problem. I have read your paper on "An Introduction to Programming the Meshless Element Free Galerkin Method," and running the codes. Would mind to help me in understanding the codes in depth? 


Adam Othman 

I'm happy to do this Adam.  Perhaps we will discover the bug in the code in the process. 

Dear Prof Dolbow,

In the coding, what is the used of connectivity array? I've run the program and have a look at the array after running it. The values are not really in inrceasing or in decreasing pattern.

One more thing, why numq2=numcell*quado^2? Do we have 16 quadrature points in each cell? But, it seems from the comment in the coding, 4 quadrature points in each cell. 



The connectivity array is used to define the background integration cells. Each column is the list of node numbers for the quadrilateral cell.   

The quadrature rule is a 2D Gauss rule that is based on two 1D rules.  So a 4-pt rule gives rise to 4x4 = 16 Gauss quadrature points.   

Dear Prof Dolbow,

Thanks for your answers. They really help my understanding of your codes. Regarding the subject above, what is the defination of tractions actually? At first I thought it it a force which means it has Newton unit. In your coding, you use tractions as shear stress

ty= -(P/(2*Imo))*((D^2)/4-gpos(2,1)^2);

which in has N/m unit. Can you explain how it is work? Because I'm thinking of applying a point load at the middle of a rectangular bar instead of load per unit length.





  Tractions have units of force per unit area.  In this example, P has units of force, I  (the area moment of inertia) has units of length^4, and D has units of length. So, if you carry everything through on the right side of the equation, you should find that the units are F / L^2.  

I also just want to correct your statement: the tractions are not used as a shear stress.  They are applied as a boundary condition.   

If you want to apply a point load, this is a little tricky with meshfree methods.  The force f_{iy} in the y-direction at each node is found from the integral over the boundary \Gamma

f_{iy} = \int_\Gamma t_y \phi_i  

where \phi_i is the nodal shape function.  So, if we think of a point load P at some point y_p on the boundary, in terms of distributions (with Dirac delta),  we obtain

f_{iy} = P \phi_i(y_p)

So, at the point where you want to apply the load, you have to first find all of the shape functions that are non-zero at that point.  Then you multiply each of these by the magnitude P of the point load.  


Dear Prof John !

My name is Van.I am very interested about your program.

I have a rectangle plate problem as follow :

  + The right edge is apllied with distribution load f_x

  + The top edge is apllied with distribution load f_y

  + The left edge is fixed with BC : u_x = 0

  + The bottom edge is fixed with BC : u_y = 0

Could help me about boundary condition of this problem ?

 In your program, f and q vector are integrated along one edge.But in my problem, i have to integrate them on two edge. 

How to integrate f and q vector for this problem ?



CanhLe's picture


Dear Prof. Dolbow,

I am studying Meshless methods and try to apply the methods for plastic problems. In order to understand meshless methods, I implement EFG for elasticity beam, but there is a problem with support domain. By taking radius of support domain dm = dmax*nodal spacing, I get good results when the number of nodes (nno) is not much. When nno is greater than 400 for the beam with 10 meters in length, bad results I get. Surprisingly, I keep dm constant (dm = 5) for different nno (including nno = 801 case), the solutions converge steadily.

What do you think about this situation? Is there the lower bound for dm?

Coming back your 1D code, I have found that when taking nno=100 (x = [0.0:.01:1.0]; of course gg = -.05:.01:0.95; gg(1) = 0.0;), singular problem appear (Warning: Matrix is singular to working precision in Matlab).

How can I overcome this problem? 

I looking for your and mechanicians' reply. 

Best wishes,

Canh V. Le, PhD Student
Department of Civil and Structural Engineering
The University of Sheffield
Mappin Street,
Sheffield, S1 3JD, UK.
Phone: GB +44 1142225724


First I apologize for the delay responding to this post.  I had missed it earlier for some reason.

One should not see the types of problems you describe with support sizes with EFG, assuming everything is done properly.  You might check to make sure that the shape functions are linearly precise.  In other words, they should sum to 1.0 at every point, and they should exactly represent a linear field in the x and y directions.  

Concerning the 1D code, the change you made to gg is not correct.   Use

gg = -.005:.01:0.995

vivek varshney's picture

hi sir

i am working in soil structure interaction problems. can i use mess free method in that like FEM. will please tell me for what type of problems  EFG is applicable.




vivek varshney 


Of course any problem that can be solved with finite elements can also be solved with EFG.  The question is whether or not there is a significant advantage to using EFG.

If the problem you're interested in involves large deformations, then there may be an advantage to using meshfree methods like EFG. 

vivek varshney's picture

thanks a lot for yr reply.

Dear Prof. Dolbow,

My name is Minh. 

In your EFG program, q is a vector obtained by integrating along the essential boundary. Also in your EFG program,  you used exact displacements, which you got from the analytic solution, to integrate.

However, if the exact solution have not been known, how can I compute vector q?

For example, the displacement along x direction of the essential boundary is restricted, that means, u_x = 0 along the essential boundary and u_y != 0 (not equal to 0).


Using FEM, I can just input u_x = 0 along the essential boundary and do not have to worry about u_y. But in EFG, how can I compute vector q, if I do not know the analytic solution of displacement?


N.Minh Nguyen 


The vector q comes from the constraints that are imposed by the Lagrange multiplier (read Section 5 of the paper that describes the code).  It assumes that you want to fix both u_x and u_y components of the displacement field along the left edge of the domain.  

If you want to make both of these components zero, then the change is simple: the q vector is simply the zero vector.  

However, if you want to only fix one component of the displacement field, you will have to modify the code quite a bit. This is because not only will the q vector need to change, but you will also need to change the G matrix and the number of discrete Lagrange multipliers.  

The code is not at all designed to do research: it is far too specific.  We had to keep it very short for it to fit into a paper.   It was simply meant to help people understand the method.   

Dear Prof. Dolbow,

I've been thinking of modifying your EFG code.

The discrete equations: 

[K G; G' 0] * [u; lambda] = [f; q]

If u_x is fixed and u_y is free, that means u_y should be treated as unknown variables. Hence, q (vector), lambda (vector) and G (matrix) should be modified.

My idea:

   q is now replaced by q_x, which is a zero vector (because u_x is fixed)

   the size of q, lambda, and G reduced (because u_y is not taken into account)


Is this idea correct? It seems quite easy. However,  when I modified the code applying this idea, I encountered some problems about matrix dimensions





What you propose is reasonable.  There are likely many places in the code where changes will need to be made, and it's possible you've just missed one.

I also want to comment that you will need to fix an additional degree of freedom somewhere.  Just constraining the x-component of displacement along one edge is not sufficient to prevent a rigid body translation.   

Thank you, Prof. Dolbow !


I will check my code again. In fact, my problem requires fixing u_x along one edge and u_y along another edge. But I think, I should fix u_x successfully along one edge first, then u_y will be similar.

I'll be back when I  finish this modification, or when other problems occur.



vivek varshney's picture

i want to know that 

 1)  the formulation for the 2-d problem will remain same for any problem. is it not depend upon the type of problem?

 2) i am using EFG for the soil material taking  elasto plastic model . how will i handle that problem?

3) in my problem there are 4 interfaces. how will i take care of interface. i mean to say at the interface how should i take nodes.And how should i assemble my matrix? how will i formulate my problem?

4) is there any mfree package for that to analysis for elasto plastic model for soil?  

kindly clear my doubt 






If you're seeking to modify the EFG Matlab routines to address your problem, I don't think it is a good idea.  As I've mentioned here before, the code is not meant for research.  It is designed for a very specific problem while being clear enough to allow people to learn the method.

Of course the formulation depends on the type of problem and governing equations you are solving.  This is the case for any numerical method. 

EFG for elasto-plasticity is really not any different than FEM for elasto-plasticity.  

Your third question is very general and I'm not sure this forum is the best place to address it.  It might be best for you to sit down with your advisor and work through these issues.

I'm not aware of any freely-available meshfree package that contains a soil model.  

Anand V Kulkarni's picture

Prof. Dolbow

Thanks a lot. I have  a few more doubts.


Is  value of shape function at point 1 in the DOI of point 2 equal to  value of the shape function at point 2 in the DOI of point1?


How does EFG method perform in the solution of non-linear problems?







I don't think one can expect symmetry of the shape functions as you've asked. It may happen for some very specific cases (uniform grids and support multipliers), but in general it is not something that is enforced and I don't see why one should expect it.

EFG is a Galerkin method (thus the 'G' in the name).  It can be applied to nonlinear problems just like the finite element method.  It can handle larger deformations than the FEM, but this comes at a greater computational cost.   

vivek varshney's picture

prof. dolbow

please give suggestion to me for queries which i have asked earlier



vivek varshney 

vivek varshney's picture

thanks for your reply

Dear Prof Dolbow !

I am using EFG for elasto plastic problems.I have a problem :

When i use Newton Raphson method for EFG i have equations

[K G;G' 0]*{du;dlamda}={f_ext - f_int ; q_ext - q_int} .

Is it right ? How could i get "f_int" and "q_int" ? Are they similar the process  of "f" and "q" in your program ?

Regards !

Long Van 

vivek varshney's picture

i want to know how to treat point load on beam in EFG.earlier you have expained. but i could not get you. please explain

You need to find all active (i.e. non-zero) EFG shape functions at the location where the point load is applied. For the sake of argument, let's assume that there are five non-zero shape functions at the point in question, and their values are [0.1 0.2 0.4 0.2 0.1].  Then multiply the magnitude of the point load by the magnitude of each of the shape functions, so this gives a vector F*[0.1 0.2 0.4 0.2 0.1].  The entries to this vector then need to be placed in the appropriate slots in the global force vector.  

vivek varshney's picture

thanks a lot sir for your reply



Dear Prof.Dolbow

I have some questions about G matrix,qk and f vector.

1/.In order to form f vector in your program,we must loop on Gauss points on traction boundary.The results are values fx,fy at nodes on traction boundary.But I think because of DOI of Gauss points so there are some values fx,fy at nodes which are not on traction boundary in f vector.This contrast with FEM,just only nodes on traction have fx,fy and others are equal zeros.

2/.G matrix and qk vector in your program are formed by looping on 1 displacement boundary.If there are 2 or 3 boundaries which are not continuous,what should I do ? How could I treat nodes directly without looping on Gauss points ?

3/.In case essential BCs are nonlinear or general,is the discrete equations still hold ?

Please explain these things to me !

Regards !


Dear Colleagues,

My name is Metsis Panagiotis and I am a doctoral candidate in the school of Civil Engineering of the National Technical University of Athens (Greece). I am trying to fully understand some of the meshless techniques and especially the Element Free Galerkin Method. For that reason I am developing a computer code (2D for now), trying to solve some simple (for the time being) problems. I was wondering if someone could answer (some of) my questions:

1) Incorporation of nodal constraints. When I am trying to solve a model with nodal constraints I come up with a very sparse G matrix (with the majority of the elements being zeros) that is problematic to invert (as part of the m matrix). Finally I acquire a solution, but with large errors (over the displacements field first of all). What should the overall approach be?

2) Background mesh. Are there any other integration methods for the EFG rather than constructing a background mesh of structured cells?

3) Shear locking. In specific models I have encountered something similar to shear locking phenomena. Is this a reality, or should I check my code again?

4) I acquire much different results with small change of variables (such as background cells or nodes). Is this true for the method?

5) My future goal is to simulate some fracture phenomena using the general pattern of meshless methods. Which of the existing meshless methods do you consider to be the most appropriate?

Thank you very much for your time.

Metsis Panagiotis

Civil Engineer, MSc

Doctoral Candidate

School of Civil Engineering

National Technical University of Athens


Alejandro Ortiz-Bernardin's picture

Dear Panagoitis,

For your sparse matrix try TAUCS. I have
successfully used it to manage large sparse matrices. It
is easy to use but has one problem. There are some OS that really
don't like TAUCS and so problems may arise. One of those systems is
Windows. But if you are familiar with UNIX/LINUX there shouldn't be
any problem. I use Red Hat Enterprise 9 and it runs pretty well.
TAUCS is open source code written in C. If you are using C++ you will
have to make a wrapper in order to compile it. If you use C, no wrap
is needed. Some of functions in Matlab and Mathematica were
implemented using TAUCS, so that should be a signal that TAUCS is not
a bad option.

Here is the link


Alejandro A.


abshaw's picture

Dear Panagiotis,

let me try to answer few of your questions.

2) nodal integration technique does not require any background mesh. you may refer to this paper  (A stabilized conforming nodal integration for Galerkin mesh-free methods, IJNME, 50, 435 - 466, 2001) and the references therein.

3) EFG also suffer  shear locking unless explicitly treated.

4) if you do integration using background mesh then slight change in mesh may affect the solution. the background mesh is required to be chosen such that the support of the kernel and background cell are aligned perfectly.  you may find some details in this paper (Numerical integration of the Galerkin weak form in meshfree methods,  Computational Mechanics,23,219 -230,1999  ).

5) I find this (Cracking particles: a simplified meshfree method for arbitrary evolving cracks, IJNME< 61, 2316 - 2343, 2004) paper very promising. 

1) i am not sure. you may try  different solver which handles sparse matrix.


this response is not comprehensive but hope it helps.



Thank you very much for your responses. The sparse matrix m (GG) I come up with, is propably due to the wrong implementation of the nodal constraints. I use the Langrange multiplier method but I can't find any analytical papers about that specific (nodal constraints) issue. Any ideas?

Thanks in advance,


Alejandro Ortiz-Bernardin's picture

have you looked at the original paper of EFG?

Alejandro A. Ortiz

This was one of the first papers I red. From my last post, I solved some of the problems mentioned there and I am expanding my code... : )

Metsis Panagiotis

gohelvivek's picture

Hello Sir,

You have been very helpful to me thankyou for your support.

Sir, can you please tell me where can i found the derivation of the exact solution which is given by Dr. Dolbow for 1d EFG problem.



In your EFG program, q is a vector obtained by integrating along the essential boundary. Also in your EFG program,  you used exact displacements, which you got from the analytic solution, to integrate.why you used the exact solution?!!!!!!!!!

when i have exact solution in my problem ,why use meshless method?

please response my problem

When setting up a problem to look at convergence studies, you have to set the bcs to match the exact solution.  Of course, one does not need to know the exact solution to use EFG (or any other numerical method) in general.

I can't understand.Prof.Dolbow used exact solution in q.If he didn't know  exact solution,how he could write q? ?

I look forward for your answers! 

Regards !

Ali Ramesh 


  In elastostatics problems, the boundary consists of two parts.  On one part of the boundary, we know the traction field.  On the other part of the boundary, we know the displacement field.  

  All you need to know to construct the vector q is the displacement field on the displacement boundary and the MLS shape functions there.  You do not need to know the exact solution to the system of equations.  You only need to know the boundary conditions, the same requirement for any boundary value problem.

I want program thin plate for natural frequencies with EFG with penelty method.I have some questions:

1-Is four gauss point enough for programming?

2-I want to find natural frequencies for Fully free thin plate.(ffff)I think for this boundary conditions i dont need use penalty factor.Is it true?(because we dont have to impose b.c there fore we dont need penalty factor)

3-I can not write program for finding gauss point for 10 gauss point in 2D problems.but i saw a program for 4 gauss can i write program for more than 4 gauss point?

Has any one program for thin plate?(with EFG or other methods)

(excuse me for my writing!)

Thanks in advance,


Hello all,

thanks for your response! I solved my problems myself.

gohelvivek's picture

Hello Dr. Dolbow,

The paper has a very nice descrition of the EFG method it is very helpful to understanding it.

How did you find the exact solution of the 1d bar problem. Please recommend if any letrature is there.


Vivek Gohel

PG student

Nirma university, India.


I think I simply derived it analytically.

 -Prof. Dolbow

Dear Prof. Dolbow,

                            i am raghavendra masters student from india . sir i went through your EFG 1D code, and i am trying to apply it to 1d heat transfer problem with internal heat generation .i am trying to reproduce paper in which element free galerkin methood is applied to heat transfer.

                             sir  my problemstatement is a bar which has 373k at one end and 200k at other end and it has uniform internal heat generation

sir i have made following changes to code to solve problem

 k1=10; area=1.0;% using k1 conductivity instead of E yougs modulus

qh=100;% factor for internal heat generation

k = k+(weight*k1*area*jac)*(dphi'*dphi)   % STIFF NESS matrix

   fbody = area*xg*qh

   f = f+(weight*fbody*jac*qh)*phi' % force matrix

instead of E i used K in stiffnes matrix formulation

in force formulation i have used qh internal heat generation 



q = [373] 


 i used 373k as temperature at one end and i want to apply 2ook at the other end

i want to know how to apply 200k at other end of bar ?

 i understood whole code but i did not get how to apply boundary condition through lagrange multipliers ,sir can you please  explain application of boundary condition

sir my mail id is please send me the paper 

 J. Dolbow and T. Belytschko (1998), "An Introduction to Programming the Meshless Element Free Galerkin Method," Archives of Computational Methods in Engineering, vol. 5, no. 3, pp. 207--242.

following are my queries

1) are modifications to code correct,when i ran code i did not get any error?

2) i want to apply 1 more boundary condition i.e. 200k at other end of bar

3) more explanation on enforcement of boundary condition by lagrange multipliers 



thanx in advance

b raghavendra

 Masters student india



satish diwakar's picture

Respected sir,

                   I am Satish diwakar PG student , my question is how to change no of nodes in your 2 d matlab code which is written for TIMESHENKO beam. i have studied your "A MESH-FREE METHOD FOR STATIC AND FREE VIBRATION ANALYSES OF THIN PLATES OF COMPLICATED SHAPE " paper . sir if you have matlab code for this paper please send on my mail ID as If any body know this kindly reply or  send mail.

                                        "Thanks in advance"

Dear Dr. J. Dolbow


Have you any  MATLAB subroutin for 3Dnode generation? i cant achive connectivity of nodes correctly.

can you help me?


Subscribe to Comments for "EFG Matlab Routines"

More comments


Subscribe to Syndicate