Home

Partial Differential Equation Toolbox User's Guide

1. c 1 0 0 0 0 0 0 c 1 0 0 0 0 0 0 c 1 0 0 0 0 0 c 1 0 0 0 0 c 1 0 0 0 0 c 1 c 1 0 0 0 0 c 2 0 0 0 0 c 1 0 0 0 0 c 2 0 0 0 0 SET 0 0 0 0 c 2 assempde The bullet symbol below means that the matrix is symmetric nc 3 c 1 c 2 0 0 0 0 c 3 0 0 0 0 0 0 c 1 c 2 0 0 0 0 c 3 0 0 0 0 0 0 c 1 c 2 0 0 0 0 c 3 n 4 c 1 c 3 0 0 0 0 c 2 c 4 0 0 0 0 0 0 c 1 c 3 0 0 0 0 c 2 c 4 0 0 0 0 0 0 c 1 c 3 0 0 0 0 c 2 c 4 The casen 3 takes precedence over the case n N and the form n N thus cannot be used n 6 c 1 0 0 0 0 0 0 c 2 0 0 0 0 0 0 c 3 0 0 0 0 0 0 c 4 0 0 0 0 0 0 c 5 0 0 0 0 0 0 c 6 5 25 assempde 12 18 20 21 10 16 19 0 c 10 5 26 assempde ne 36 c 1 c 3 13 15 25 27 c 2 c 4 14 16 26 c 28 c 5 c 7 17 19 29 c 31 c 6 c 8 18 20 30 c 32 c 9 c 11 21 23 33 c 35 c 10 c 12 22 24 34 c 36 See Also initmesh refinemesh pdebound assema assemb 5 27 csgchk Purpose Synopsis Description See Also 5 28 Check validity of Geometry Description matrix gstat csgchk gd xlim ylim gstat csgchk gd gstat csgchk gd xlim ylim checks if the solid objects in the Geometry Descriptio
2. Examples of Elliptic Problems The problem can also be solved by typing gt Compare with solution not using subdomains K F assempde Ishapeb p e t 1 0 1 u1 K1F norm u ul inf pdesurf p t u Y ou can run this entire example by typing pdedemo4 0 15 0 1 0 05 i ay iy ii IN F a a A N EQ AIN O NN ES mh ANUN NNA 0 5 pa He Wy MN ANN oe U a UD pd o 2 15 2 Examples Examples of Parabolic Problems This section describes the solution of some parabolic PDE problems The problems are solved using both the graphical user interface and the command line functions of the PDE Toolbox The Heat Equation A Heated Metal Block A common parabolic problem is the heat equation du d5 Au 0 The heat equation describes the diffusion of heat in a body of some kind See the section Application M odes for more information about heat transfer and diffusion problems This first example studies a heated metal block with a rectangular crack or cavity The left hand side of the block is heated to 100 degrees centigrade At the right hand side of the metal block heat is flowing from the block to the surrounding air at a constant rate All the other block boundaries are isolated This leads to the following set of boundary conditions when proper scaling of t is chosen e y 100 on the left hand side Dirichlet condition o 10 on the right
3. Basis 10 Time 2 70 ew conv eig 0 Basis 13 Time 3 50 ew conv eig 0 Basis 16 Time 4 36 ew conv eig 0 Basis 19 Time 5 34 ew conv eig 1 Basis 22 Time 6 46 ew conv eig 2 Basis 25 Time 7 61 ew conv eig 3 Basis 28 Time 8 86 ew conv eig 3 Basis 31 Time 10 23 ew conv eig 5 Basis 34 Time 11 69 ew conv eig 5 Basis 37 Time 13 28 ew conv eig 7 Basis 40 Time 14 97 ew conv eig 8 Basis 43 Time 16 77 ew conv eig 9 Basis 46 Time 18 70 ew conv eig ll Basis 49 Time 20 73 ew conv eig ll Basis 52 Time 22 90 ew conv eig 13 Basis 55 Time 25 13 ew conv eig 14 Basis 58 Time 27 58 ew conv eig 14 Basis 61 Time 30 13 ew conv eig 15 Basis 64 Time 32 83 ew conv eig 16 Basis 67 Time 35 64 ew conv eig 18 Basis 70 Time 38 62 ew conv eig 22 End of sweep Basis 70 Time 38 62 ew conv eig 22 Basis 32 Time 43 29 ew conv eig 0 Basis 35 Time 44 70 ew conv eig 0 Basis 38 Time 46 22 ew conv eig 0 Basis 41 Time 47 81 ew conv eig 0 Basis 44 Time 49 52 ew conv eig 0 Basis 47 Time 51 35 ew conv eig 0 Basis 50 Time 53 27 ew conv eig 0 Basis 53 Time 55 30 ew conv eig 0 End of sweep Basis 53 Time 55 30 ew conv eig 0 2 30 Examples of Eigenvalue Problems You can see that two Arnoldi runs were made In the first 22 eigenvalues converged after a basis of size 70 was computed in the second where the vectors were
4. a ST Save As Filter home pdedemo m Directories Files home pdedemo es 3 home pdedemo ar Lia A Selection home pdedemoj Done Filter Cancel _ 5 Save As displays a dialog box in which you can specify the name of the file in which to save the CSG model and other information regarding the GUI session Y ou can also change the directory in which it is saved If the filename is given without a m extension m is appended automatically The GUI session is stored in a Model M file which contains a sequence of drawing commands and commands torecreate the modeling environment axes scaling grid etc If you have already defined boundary conditions PDE coefficients created a triangular mesh and solved the PDE further commands to recreate the modeling and solution of the PDE problem are also included in the Model M file Thepdet oo GUI can be started from the command line by entering the name of a Model M file The model in the file is then directly loaded into the GUI 3 The Graphical User Interface 3 6 PDE Toolbox Untitled Print Device Option dps Paper Orientation portrait 4 Send to Printer Y Send to File Help Cancel Print al Print displays a dialog box for printing a hardcopy of the PDE Toolbox figure Only the main part of the figure is printed not the upper and lower menu and information parts In
5. 000 cece cece ee eeaee 3 3 EME MU A A A a ia 3 3 EOI MUI a ee 8 Sean 3 7 Options Meli e ea 3 9 Draw Menu ocococcccocc eee nett eens 3 13 ROtHtGs ia bee bea HERS Oe oe ea aaa 3 14 Boundary Menu 00 eee eee 3 15 PDE Men wis Saker A ie ae Pe ee et 3 18 MieshiMeAU te a dd bad 3 22 SOIVEMENU i erri nie kea eee eee eee eens 3 25 POE MCA Lati Sok Groton a ta a howto Rated 3 30 Window Menu drae naa ea Ra eee eee ees 3 37 Help Mend e Me blew thd Alana gE see 3 37 The TOOIDAF 0000 o ci A a 3 38 3 The Graphical User Interface ThePDE Toolbox includes a graphical user interface GUI pdetoo l The main components of the GUI arethe menus the dialog boxes and the toolbar They are explained in detail in this section PDE Toolbox M enus PDE Toolbox Menus The GUI has a pull down menu bar that you can use to control the modeling It conforms to common pull down menu standards Menu items followed by a right arrow lead to a submenu Menu items followed by an ellipsis lead toa dialog box Stand alone menu items lead to direct action Some menu items can be executed by using keyboard accelerators pdetool also contains a toolbar with icon buttons for quick and easy access to some of the most important pdetoo functions The following sections describe the contents of the pdet oo menus and the dialog boxes associated with menu items File Menu File Edit Options D New lt Ctrl gt n Open lt Ctrl gt
6. 4 29 4 The Finite Element M ethod Solving Tw f using the discrete sine transform would not be an advantage in itself since the system is tridiagonal and should be solved as such However for the full system Kv F a transformation of the blocks in K turns it into N 1 decoupled tridiagonal systems of size N 1 Thus a solution algorithm would look like 1 DivideF into N 1 blocks of length N 1 and perform an inverse discrete sine transform on each block 2 Reorder the elements and solve N 1 tridiagonal systems of size N 1 with 2 h h h h 2h h cos xi N on the diagonal and h h on the subdiagonals 3 Reverse the reordering and perform N 1 discrete sine transforms on the blocks of length N 1 When using a fast solver such as this one time and memory are also saved since the matrix K in fact never has to be assembled A drawback is that since the mesh has to be regular it is impossible to do adaptive mesh refinement The fast elliptic solver for Poisson s equation is implemented in poi solv The discrete sine transform and the inverse discrete sine transform are computed by dst andi dst respectively 4 30 Reference Commands Grouped by Function 5 3 PDEAIJO It MS osn iar ths Seo a yale er WR Rae a 5 3 User InterfaceAlgorithms 00 cee eee 5 3 Geometry Algorithms 0 00 ee 5 4 Plot FUNCtONS 23 va wel e Ede oe cow Eee ee oe Se eee 5 4
7. The braced quantity is the jump in normal derivative of v across edge q h is the length of edge 1 and the sum runs over Ej the set of all interior edges of the triangulation The coefficients a and are independent of the triangulation This bound is turned into an element wise error indicator function E K for element K by summing the contributions from its edges The final form for the toolbox equation V cVu au f becomes 1 2 EK ah awl B A y htm cu 1 te OK where n is the unit normal of edge t and the braced term is the jump in flux across the element edge The L norm is computed over the element K This error indicator is computed by the pdej mps function The Mesh Refiner The PDE Toolbox is geared to elliptic problems For reasons of accuracy and ill conditioning they require the elements not to deviate too much from being equilateral Thus even at essentially one dimensional solution features such as boundary layers the refinement technique must guarantee reasonably shaped triangles 4 27 4 The Finite Element M ethod When an element is refined new nodes appear on its midsides and if the neighbor triangle is not refined in a similar way it is said to have hanging nodes The final triangulation must have no hanging nodes and they are removed by splitting neighbor triangles To avoid further deterioration of triangle quality in successive generations the longest edge bisection scheme Rose
8. p e t refinemesh Ishapeg p e t v pdeeig Ilshapeb p e t 1 0 1 Inf 100 1 first eigenvalue pdesurf p t v 1 first eigenmode figure membrane 1 20 9 9 the MATLAB function figure 1 16 sixteenth eigenvalue pdesurf p t v 16 sixteenth eigenmode In the standard case c and d are positive in the entire region All eigenvalues are positive and 0 is a good choice for a lower bound of the interval The cases where either c or d is zero are discussed below e If d 0 in a subregion the mass matrix M becomes singular This does not cause any trouble provided that c gt 0 everywhere The pencil K M has a set of infinite eigenvalues e fc 0in a subregion the stiffness matrix K becomes singular and the pencil K M has many zero eigenvalues With an interval containing zero pdeei g goes on for a very long time to find all the zero eigenvalues Choose a positive lower bound away from zero but below the smallest nonzero eigenvalue e f there is a region where both c 0 and d 0 we get a singular pencil The whole eigenvalue problem is undetermined and any value is equally plausible as an eigenvalue Some of the awkward cases are detected by pdeei g If the shifted matrix is singular another shift is attempted If the matrix with the new shift is still singular a good guess is that the entire pencil K M is singular If you try any problem not belonging to the standard
9. y c Vu V c Vu yy yoy fy y c21 Vu V Ca Vu 47 47 azzi h You can work with systems of arbitrary dimension from the command line F or the elliptic problem an adaptive mesh refinement algorithm is implemented It can also be used in conjunction with the nonlinear solver In addition a fast solver for Poisson s equation on a rectangular grid is available 1 Tutorial The following boundary conditions are defined for scalar u e Dirichlet hu r on the boundary Q e Generalized Neumann cVu qu g OndQ is the outward unit normal g q h and r are complex valued functions defined on o The eigenvalue problem is a homogeneous problem i e g 0 r 0 In the nonlinear case the coefficients g q h and r can depend on u and for the hyperbolic and parabolic PDE the coefficients can depend on time For the two dimensional system case Dirichlet boundary condition is hjju hy74 r haju hoi 1 the generalized Neumann boundary condition is c11 Vu c Vu 9114 9124 8 Ivy a C21 Vuy h c7Vu 9714 9774 82 Ivy and the mixed boundary condition is hjit hizi r R Cy Vu C13 Vu Gy My qizu 8 tye Iv Xy E C21 Vu A c7 Vu 99141 994 82 hih where u is computed such that the Dirichlet boundary condition is satisfied Dirichlet boundary conditions are also called essential boundary conditio
10. 0 5 x AN NN N SIS N N N N N N z The error function pdedemol performs all the above steps 2 Examples 2 6 A Scattering Problem The scattering problem is to compute the waves reflected from an object illuminated by incident waves F or this problem consider an infinite horizontal membrane subjected to small vertical displacements U The membrane is fixed at the object boundary We assume that the medium is homogeneous so that the wave speed is constant c NOTE Do not confuse this c with the parameter c appearing in the PDE Toolbox When the illumination is harmonicin time we can compute the field by solving a single steady problem With U x y t u x y e the wave equation turns into o 2u Au 0 or the Helmholtz s equation Au ku 0 Examples of Elliptic Problems where k the wavenumber is related to the angular frequency o the frequency f and the wavelength A by pe 2M _ 21 c c A We have yet to specify the boundary conditions Let the incident wave be a plane wave traveling in the direction cos a sin a i ka ot Pas x V x y 1 v x ye ika where v x y e u is the sum of v and the reflected waver u v r The boundary condition for the object s boundary is easy u 0 i e r v x y For acoustic waves where v is the pressure disturbance the proper condition would be ou 0 Thereflected waver travels outward from the obje
11. 2 Examples 2 38 and u in turn is defined by v l v u 2G are volume forces Thisisan elliptic PDE of system type u is two dimensional but you need only toselect the application mode Structural Mechanics Plane Stress and then enter the material dependent parameters E and v and the volume forces k into the PDE Specification dialog box In this mode you can also solve the eigenvalue problem which is described by V c Vu Adu mat p the density can also be entered using the PDE Specification dialog box In the Plot Selection dialog box the x and y displacements u and v and the absolute value of the displacement vector u v can be visualized using color contour lines or z height and thedisplacement vector field u v can be plotted using arrows or a deformed mesh In addition for visualization using color contour lines or height you can choose from 15 scalar tensor expressions Ox _ du e ux exx the x direction strain ex eyy the y direction strain ey Application M odes exy the shear strain Yxy sxx the x direction stress ox syy the y direction stress oy sxy the shear stress tyy el the first principal strain e1 e2 the second principal strain e gt 51 the first principal stress 01 s 2 the second principal stress o gt e von Mises the von Mises effective stress os 05 0 0 For a more detailed discussion on the theory of stress strain rela
12. 3 34 option is turned off the solution plot is a 2 D plot and is plotted in the main axes ofthepdetoo GUI If the Height 3 D plot option is used the solution plot is a 3 D plot in a separate figure window If possible the 3 D plot uses an existing figure window If you would liketo plot in a new figure window simply typef igure at the MATLAB command line Additional Plot Control Options In the middle of the dialog box are a number of additional plot control options Plot in x y grid f you check this option the solution is converted from the original triangular grid toa rectangular x y grid This is especially useful for animations since it speeds up the process of recording the movie frames significantly Show mesh In the surface plots the mesh is plotted using black color if you check this option By default the mesh is hidden Contour plot levels For contour plots the number of level curves e g 15 or 20 can be entered Alternatively a MATLAB vector of levels can be entered The curves of the contour plot are then drawn at those levels The default is 20 contour level curves Examples 0 100 1000 logspace 1 1 30 Colormap Using the Colormap pop up menu you can select from a number of different colormaps coo gray bone pink copper hot jet hsv and prism Plot solution automatically This option is normally selected If turned off the PDE Toolbox does not display a plot of the solution immediately upon solvi
13. Setformula co de EOS EE Ba The adaptively refined mesh The solution of the AC power electromagnetics equation is complex The plots show the real part of the solution a warning message is issued but the solution vector which can be exported to the main workspace is the full complex solution Also you can plot various properties of the complex solution by using the user entry i mag u andabs u aretwo examples of valid user entries The skin effect is an AC phenomenon Decreasing the frequency of the alternating current results in a decrease of the skin effect Approaching DC conditions the current density is close to uniform experiment using different angular frequencies Color i Height i The current density in an AC wire Application Modes Conductive Media DC For electrolysis and computation of resistances of grounding plates we have a conductive medium with conductivity o and a steady current The current density J is related to the electric field E through J oE Combining the continuity equation V J Q where Q is a current source with the definition of the electric potential V yields the elliptic Poisson s equation V oVV Q The only twoPDE parameters are the conductivity o and the current source Q The Dirichlet boundary condition assigns values of the electric potential V to the boundaries usually metallic conductors The
14. ag jo dx f atas WY Jal foie f s jas i Ap Ny 4 5 4 The Finite Element M ethod Use the following notations Kegs Jv Vo dx Stiffness matrix M j aide Mass matrix Qij f 10 0 ds F foie G fastas and rewrite the system in the form K M Q U F G K M and Q are N p by Np matrices and F and G are N pvectors K M and F are produced by as sema while Q G are produced by as semb When it is not necessary to distinguish K M and Q or F and G we collapse the notations to KU F which form the output ofassempde When the problem is sef adjoint and el lipticin the usual mathematical sense the matrix K M Q becomes symmetric and positive definite M any common problems have these characteristics most notably those that can also be formulated as minimization problems For the case of a scalar equation K M and Q are obviously symmetric If c x gt 8 gt 0 a x gt O and q x gt O with q x gt 0 on some part of dQ then if U FO U K M Q U c Vul au dx qu ds gt 0 if U 0 Q dQ UT K M Q U isthe energy norm There are many choices of the test function spaces The toolbox uses continuous functions that are linear on each triangle of the mesh Piecewise linearity guarantees that the integrals defining the stiffness matrix K exist Projection onto Vy is nothing more than linear interpolation and the evaluation of the solution inside a triangle is done just in terms
15. try with a smaller interval spd is 1 if the pencil is known to be symmetric positive definite default 0 tol conv isthe expected relative accuracy Default is100 eps whereeps isthe machine precision j max is the maximum number of basis vectors The algorithm needs j ma x n working space soa small value may be justified on a small computer otherwise let it be the default value max 100 Normally thealgorithm stops earlier when enough eigenvalues have converged max mul is the number of Arnoldi runs tried Must at least be as large as maximum multiplicity of any eigenvalue If a small value of j max is given many Arnoldi runs are necessary The default value is max mu n which is needed when all the eigenvalues of the unit matrix are sought sptarn Algorithm Cautionary TheArnoldi algorithm with spectral transformation is used The shift is chosen atub l b or at arandom point in interval b ub when both bounds are finite The number of stepsj in the Arnoldi run depends on how many eigenvalues there are in the interval but it stops at j mi n j max n After a stop the algorithm restarts to find more Schur vectors in orthogonal complement to all those already found When no more eigenvalues arefoundin b lt mb lt ub the algorithm stops For small values of j max several restarts may be needed before a certain eigenvalue has converged The algorithm works when j max is at least one larger than the number of
16. ud atan cos pi 2 x gt ut0 3 sin pi x exp sin pi 2 y n 31 tlist linspace 0 5 n list of times Y ou are now ready to solve the wave equation The general form for the hyperbolic PDE in the PDE Toolbox is y d V cVu au i ot so here you haved 1 c 1 a 0 and f 0 uu hyperbolic ud ut0 tlist squareb3 p e t 1 0 0 1 2 25 2 Examples To visualize the solution you can animate it Interpolate toa rectangular grid to speed up the plotting delta 1 0 1 1 uxy tn a2 a3 tri2grid p t uu 1 delta delta gt gp tn a2 a3 umax max max uu umi n mi n min uu gt newpl ot M moviein n for i l n pdeplot p e t xydata uu i zdata uu i mesh off xygrid on gridparam gp colorbar off zstyle continuous axis 1 1 1 1 umin umax caxis umin umax M 1 getframe end gt movie M 10 Y Y ou can find a complete demo of this problem including animation in pdedemo6 If you have lots of memory you can try increasing n the number of frames in the movie Animation of the solution to the wave equation 2 26 Examples of Eigenvalue Problems Examples of Eigenvalue Problems This section describes the solution of some eigenvalue PDE problems The problems are solved using the graphical user interface and the command line functions of the PDE Toolbox Eigenvalues and Eigenfunctions
17. o Save As Print Exit lt Ctrl gt w New Create a new empty Constructive Solid Geometry CSG model Open Load a model M file from disk Save Save the GUI session to a Mode M file Save As Save the GUI session to a new Model M file Print Print a hardcopy of the PDE Toolbox figure Exit Exit thepdet ool graphical user interface New New deletes the current CSG model and creates a new empty model called Untitled 3 The Graphical User Interface 34 Open Open Filter home pdedemos Directories Files 3 mm shome pdedemoy 4 lt a P Ena 4 Selection home pdedemostest m Done Filter Cancel Ftd Open displays a dialog box with a list of existing M files from which you can select the filethat you want toload You can list the contents of a different directory by changing the path in the Selection text box You can use the scrollbar to display more filenames You can select a file by double clicking on the filename or by clicking on the filename and then clicking on the Done button When you select a file the CSG model that is stored in the M odel M file is loaded intothe PDE Toolbox and displayed Also the equation the boundary conditions and information about the mesh and the solution are loaded if present and the modeling and solution process continues tothe same status as when you saved the file PDE Toolbox M enus Save As
18. xydata xystyle contour zdata zstyle flowdata flowstyle col or map xygrid gridparam mes h col orbar title levels data off flat interp off on data off onti nuous discontinuous data off rrow colormap cool pff lon tn a2 a3 pff on off on 10 Triangle data x y data plot style Show contours Node or triangle data 3 D height plotstyle Node or triangle data Flow plot style x y data colormap name or colormap matrix Convert to x y grid before plotting Triangle index and interpolation parameters from earlier call tot ri 2grid Show mesh in plot Show colorbar Plot title text Number of levels or a vector specifying levels 5 69 pdeplot Examples See Also 5 70 Thepdepl ot is used both from inside thepdet 00 GUI and from the command line It is able to display three entities simultaneously xydat a can be visualized by a surface plot Either flat or interpolated default shading can be used for the surface plots A contour plot can be superimposed on the surface plot in black or plotted independently in colors by settingcontour toon zdata is visualized by displaying height The triangles can be either tilted by interpolation default or flat Flow data can be visualized by plotting arrows likethe MATLAB qui ver plot All data types can be either node data or triangle data flow data can only be triangle data Node data is represented by a column
19. 0 with Dirichlet boundary conditions u e Y_ Tol fixes the exit criterion from the Gauss N ewton iteration i e the iterations are terminated when the residual norm is less then To The norm in which the residual is computed is selected through Norm This can be any admissible MATLAB vector norm or energy for the energy norm 5 67 pdenonlin Diagnostics See Also 5 68 Maxlter andMi nStep are safeguards against infinite Gauss N ewton loops and they bound the number of iterations and the step size used in each iteration Setting Report toon forces printing of convergence information If the Newton iteration does not converge the error messageToo many iterations orStepsize too small is displayed If theinitial guess produces matrices containing NaN or nf elements the error message Unsuitable initial guess UO default U0 0 is printed assempde pdebound pdeplot Purpose Synopsis Description Generic PDE Toolbox plot function pdeplot p e t PropertyName PropertyVal ue h pdeplot p e t PropertyName PropertyVal ue pdeplot p e t pl vl is the generic PDE Toolbox plot function It can display several functions of a PDE solution at the same time The geometry of the PDE problem is given by the mesh datap e andt Details on the mesh data representation can be found in the entry on i nit mesh Valid property value pairs include Property Name Property Value Default Description
20. 000 00 e eae 4 21 Adaptive Mesh Refinement 0 0048 4 26 The Error Indicator Function 0020 0c eee eee 4 26 The Mesh Refiner 0 00 00 cee eee 4 27 The Termination Criteria 0 00 eee eee 4 28 Fast Solution of Poisson s Equation 4 29 4 The Finite Element M ethod 4 2 The core of the PDE Toolbox is a PDE solver that uses the Finite Element Method FEM for problems defined on bounded domains in the plane The toolbox provides a user friendly interface between the computing environment of MATLAB and the technical procedures of the Finite Element M ethod The Elliptic Equation The Elliptic Equation The basic elliptic equation handled by the tool box is V cVu au f in Q where Q is a bounded domain in the plane c a f and the unknown solution u are complex functions defined on Q c can also be a 2 by 2 matrix function on Q The boundary conditions specify a combination of u and ts normal derivative on the boundary e Dirichlet hu r on the boundary 0Q e Generalized Neumann cVu qu gon aa e Mixed Only applicable to systems A combination of Dirichlet and generalized Neumann is the outward unit normal g q h andr are functions defined on dQ Our nomenclature deviates slightly from the tradition for potential theory where a Neumann condition usually refers to the case q 0 and our Neumann would be called a mixed condition In some c
21. FEM formulation of the scalar PDE problem ty cVu au f on Q or the system PDE problem Ju du Y cC Vu tau f on Q on a mesh described by p e and t with boundary conditions given by b and with initial valueu0 For the scalar case each row in the solution matrix u1 is the solution at the coordinates given by the corresponding column inp Each column in u1 is the solution at the time given by the corresponding item int ist For a system of dimension N with np node points the first np rows of ul describe the first component of u the following np rows of u1 describe the second component of u and so on Thus the components of u are placed in the vector u as N blocks of node point rows b describes the boundary conditions of the PDE problem b can be either a Boundary Condition matrix or the name of a Boundary M file The boundary conditions can depend ont the time The formats of the Boundary Condition matrix and Boundary M file are described in the entries on assemb and pdebound respectively The geometry of the PDE problem is given by the mesh datap e andt Details on the mesh data representation can be found in the entry oni nit mesh The coefficientsc a d andf of the PDE problem can be given in a variety of ways The coefficients can depend ont the time A complete listing of all options can be found in under the entry ass empde parabolic Examples See Also atol andrtol are absolute and relativ
22. K F assempde b p e t c a f assembles the PDE problem by approximating the Dirichlet boundary condition with stiff springs see the section The Elliptic System in The Finite Element Method chapter for details K andF arethestiffness matrix and right hand side respectively The solution tothe FEM formulation of the PDE problemisu K F K F B ud assempde b p e t c a f assembles the PDE problem by eliminating the Dirichlet boundary conditions from the system of linear equations u1 K F returns the solution on the non Dirichlet points The solution to the full PDE problem can be obtained as the MATLAB expression u B ul ud K M F Q G H R assempde b p e t c a f gives a split representation of the PDE problem u assempde K M F Q G H R collapses the split representation into the single matrix vector form and then solves the PDE problem by eliminating the Dirichlet boundary conditions from the system of linear equations K1 F1 assempde K M F Q G H R collapses thesplit representation into the single matrix vector form by fixing the Dirichlet boundary condition with large spring constants K1 F1 B ud assempde K M F Q G H R collapses the split representation into the single matrix vector form by eliminating the Dirichlet boundary conditions from the system of linear equations 5 19 assempde 5 20 b describes the boundary conditions of the PDE problem b can be either a Boundary Condition ma
23. Utility Algorithms esete apiye iia a aE ees 5 5 User Defined Algorithms oococcooocococrr 5 7 5 Reference This chapter contains detailed descriptions of most of the functions in the Partial Differential E quation PDE Toolbox Some of the functions that are used by the graphical user interface GUI function pdetool are not listed in this chapter because they are not supposed to be used directly they are called frompdetool 5 2 Commands G rouped by Function Commands Grouped by Function PDE Algorithms Function Purpose adaptmesh Adaptive mesh generation and PDE solution assema Assemble area integral contributions assemb Assemble boundary condition contributions assempde Assemble the stiffness matrix and right hand side of a PDE problem hyperbolic Solve hyperbolic PDE problem parabolic Solve parabolic PDE problem pdeeig Solve eigenvalue PDE problem pdenonlin Solve nonlinear PDE problem poisolv Fast solution of Poisson s equation on a rectangular grid User Interface Algorithms Function Purpose pdecirc Draw circle pdeellip Draw ellipse pdemdl cv Convert PDE Toolbox 1 0 Model M files pdepol y Draw polygon pderect Draw rectangle pdetool PDE Toolbox graphical user interface GUI 5 3 5 Reference Geometry Algorithms Function Purpose csgchk csgdel decsg initmesh jigglemesh pdearcl poi mesh refinemesh wbound wgeom Check validity of
24. assem BE E SS K M mn assem alp Dp p p t p p f f u0 F u0 ti me f u0 time sdl f time f time sdl c a f assembles the stiffness matrix K the mass matrix m and the right hand side vector F Theinput parametersp t c a f u0 time andsdl havethesame meaning as inassempde assempde 5 13 assemb Purpose Synopsis Description 5 14 Assemble boundary condition contributions Q G H R assemb b p e Q G H R assemb b p e u0 Q G H R assemb b p e u0 ti me Q G H R assemb b p e u0 ti me sdl Q G H R assemb b p e ti me Q G H R assemb b p e time sdl Q G H R assemb b p e assembles the matrices Q and H and the vectors G andR Q should be added to the system matrix and contains contributions from mixed boundary conditions G should be added to the right hand side and contains contributions from generalized Neumann and mixed boundary conditions The equation H u R represents the Dirichlet type boundary conditions The input parameters p e u0 time andsdl have the same meaning as in assempde b describes the boundary conditions of the PDE problem b can be either a Boundary Condition matrix or the name of a Boundary M file The format of the Boundary Condition matrix is described below and the format of the Boundary M file is described in the entry on pdebound The toolbox treats the follow
25. o 2 59 DIFUSIO da ea o aa tv nod Mt aE 0 RA BE AA wen 2 61 The Graphical User Interface 3 PDE Toolbox Menus 0 00 cece eee eee ees 3 3 File Mentes da aras ma Paw oP oes eee wore aa 3 3 NOW ire oho cee a a Me ie ee ee RO 3 3 OPEN ta are tis Sieur eae a eat nase os 3 4 SIMA Ss Aaa Le Ee a eS Oi ees 3 5 PRIN Gwe EEEE wae kw w haw ade oo ee med Aw Badan 3 6 EsGit MCN U a A A a as ee AO ioe Rea 3 7 Pastilla er ee ee 3 8 Options Menus erii Aa wk ee a wie 3 9 Grid Spa iNg dada e 3 10 AXeS MIES taa Bee don a dt Mae e A wel 3 11 Application taa a a AE A a eee 3 11 Draw M nu 245480438 22 w Rein O A es do 3 13 ROO oc en oe es a Wa Bie ee oa ds 3 14 Boundary Menu 0000 cece eee eee 3 15 Specify Boundary COnditi0NS o oooooooommo o 3 16 PDE M nU reran vive eres a Shea ia 3 18 PDE SpeciiGation ui 3 19 Mesh Melissa Pb a da dae 3 22 ParameterS 1 ccc eee eee eee eens 3 23 SOME MU a Sue e ie a Rie a dei 3 25 Parameter Si Ln SG BROAN a a 3 25 POEMA a a da Da ek 3 30 ParameterS 1 ccc cc ee eet eens 3 30 Additional Plot Control Options 00 3 34 Window Menu ccc eee eens 3 37 Fel DM A ranked eer ee ak ees Gee Oe 3 37 The Toolbar meso oe ii ate Bai RS EG a ES 3 38 iv Contents The Finite Element Method 4 The Elliptic Equation a 4 3 The Elliptic System 0 0 0c cee es 4 10 The Parabolic Equation 0 0c ee eeeee 4
26. refinement wherethe longest edge of each specified triangleis bisected can be demanded by giving ongest asa final parameter Usingregul ar asa final parameter results in regular refinement Some triangles outside of the specified set may also be refined in order to preserve the triangulation and its quality refinemesh Examples Algorithm See Also Refinethe mesh of the L shaped membrane several times Plot the mesh for the geometry of the L shaped membrane p e t initmesh Ishapeg hmax inf subplot 2 2 1 pdemesh p e p e t refinemesh lshapeg p e t subplot 2 2 2 pdemesh p e p e t refinemesh lshapeg p e t subplot 2 2 3 pdemesh p e p e t refinemesh lshapeg p e t subplot 2 2 4 pdemesh p e subplo The algorithm is described by the steps below 1 Pick the initial set of triangles to be refined 2 Either divide all edges of the selected triangles in half regular refinement or divide the longest edge in half longest edge refinement 3 Dividethe longest edge of any triangle that has a divided edge 4 Repeat step 3 until no further edges are divided 5 Introduce new points of all divided edges and replace all divided entries in e by two new entries 6 Form the new triangles If all three sides are divided new triangles are formed by joining the side midpoints If two sides are divided the midpoint of the longest ed
27. the button with the solution plot icon or by selecting Parameters fromthe Plot menu This opens a dialog box Boundary Condition Dialog Box In this dialog box the boundary condition for the selected boundaries is entered T he following boundary conditions can be handled e Dirichlet hu r on the boundary Generalized Neumann cVu qu g on the boundary e Mixed a combination of Dirichlet and generalized Neumann condition is the outward unit length normal The boundary conditions can be entered in a variety of ways See the entry on assemb and the Boundary Menu on page 3 15 PDE Specification Dialog Box In this dialog box the type of PDE and the PDE coefficients are entered The following types of PDEs can be handled e Elliptic PDE V cVu au f du e Parabolic PDE dV cVu au f y e Hyperbolic PDE d V cVu au f ot e Eigenvalue PDE V cVu au Adu for x and y on the problem s 2 D domain Q The PDE coefficients can be entered in a variety of ways See the entry on assempde andthe PDE Menu on page 3 18 pdetool Model M file The Modd M file contains the MATLAB commands necessary to create a CSG model It can also contain additional commands to set boundary conditions definethe PDE create the mesh solve the pde and plot the solution This type of M file can be saved and opened from the File menu The Model M file is a MATLAB function and not a script This w
28. the mid point rule This approximation is optimal since it has the same order of accuracy as the piecewise linear interpolation Consider a triangle given by the nodes P4 P2 and Pz as in the following figure P PS Xi The local triangle P1P2P3 1 The local 3 by 3 matrices contain the integrals evaluated only on the current triangle The coefficients are assumed constant on the triangle and they are evaluated only in the triangle barycenter 4 7 4 The Finite Element M ethod The simplest computations are for the local mass matrix m area AP P P3 aP JO DAA APIS M 1 gt g lagi a where P is the center of mass of A P P gt P3 i e P P gt P E E The contribution to the right hand side F is just area AP P P3 i AP 3 For the local stiffness matrix we have to evaluate the gradients of the basis functions that do not vanish on PPP Since the basis functions are linear on thetriangleP P P3 the its are constants Denote the basis functions 4 b gt and o3 sich that 0 P 1 If P P3 x y l then we have that 1 yy V ae 0 2area AP P P E and after integration taking c as a constant matrix on the triangle 1 kij 4area AP P P yer g pee of el If two vertices of the triangle lie on the boundary 0Q they contribute to the lineintegrals associated tothe boundary conditions If the two boundary points are P and P then we have OF q P yf Pal las jb i j 1 2 a
29. 0 By Uy h Ot kee S fs u f By B By C c c and the solution is given by block elimination Ela E 1_ E LE 1 1 1 C B K B BK B B3K3 B3 u f B K fi BK f2 B3K3 1 T uy K f B 4 2 13 2 Examples In the MATLAB solution below a more efficient algorithm using Choleski factorization is used time np size p 2 Find common points c pdesdp p e t nc length c C zeros nc nc FC zeros nc 1 il cl pdesdp p e t 1 icl pdesubix c cl K F assempde Ishapeb p e t 1 0 1 time 1 Kl K 11 11 d symmmd K1 1i1 i1 d Kl chol K1 d d B1 K c1 i1 a1 B1 Kl1 C ic1 ic1 C 1c1 1c1 K c1 c1 al al fl F 11 e1 K1 1f1 FC ic1 FC ic1 F c1 al el 2 c2 pdesdp p e t 2 ic2 pdesubix c c2 K F assempde Ishapeb p e t 1 0 1 time 2 K2 K 12 12 d symmmd K2 12 i2 d K2 chol K2 d d B2 K c2 12 a2 B2 K2 C ic2 i 2 C ic2 ic2 K c2 2 a2 a2 f2 F 12 e2 K2 1f2 FC 1c2 FC c2 F c2 a2 e2 13 c3 pdesdp p e t 3 ic3 pdesubix c c3 K F assempde Ishapeb p e t 1 0 1 time 3 K3 K 13 13 d symmmd K3 13 13 d K3 chol K3 d d B3 K c3 i3 a3 B3 K3 C ic3 ic3 C ic3 ic3 K c3 c3 a3 a3 f 3 F i 3 e3 K3 f3 FC ic3 FC ic3 F c3 a3 e3 Solve u zeros np 1 u c C FC u il K1 el al u cl u 12 K21 e2 a2 u c2 u i 3 K3 e3 a3 u c3
30. 1 5 How Can Solvea PDE Problem 00020 0 eee 1 6 Can Use the Toolbox for Nonstandard Problems 1 6 How Can Visualize My Results 000 e eee 1 6 Are There Any Applications Already Implemented 1 7 Can Extend the Functionality of the Toolbox 1 7 How Can I Solve 3 D Problems by 2 D Models 1 8 Getting Started 0 cece ee 1 9 Basics of The Finite Element Method 1 18 Using the Graphical User Interface 1 23 The PDE Toolbox Graphical User Interface 1 23 The Menus sci ce ave e eee ee ee ee ees 1 24 The 00lb ar unreal 1 25 The GUI Modes 0 ects 1 26 The CSG Model and the Set Formula o 1 27 Creating Rounded Corners 2 00 c eee eee 1 28 Suggested Modeling Method 00 cee eee neces 1 31 Object Selection Methods 0 0 c eee eee 1 35 Display Additional Information 000000 1 35 Entering Parameter Values as MATLAB Expressions 1 36 Using PDE Toolbox version 1 0 Model M files 1 36 Using Command Line Functions 1 37 Data Structures and Utility Functions 1 37 Hints and Suggestions for Using Command Line Function 1 40 1 Tutorial Introduction This section attempts to answer some of the questions you might formulate when you turn the first page What does this toolbox do
31. 13 The Hyperbolic Equati0M oo o ooooo 4 16 The Eigenvalue Equation ooo oooooo 4 17 Nonlinear Equations 0 0 0 0 4 21 Adaptive Mesh Refinement 4 26 The Error Indicator Function 0 0 c eee ee eee 4 26 The Mesh Refiner 0 ccc cee ene een 4 27 The Termination Criteria 0 0 asua cee eee 4 28 Fast Solution of Poisson s Equation 4 29 vi Contents Reference Commands Grouped by Function 5 3 PDE AlgorIthMS ese uscar ree eh 5 3 User Interface Algorithms ooooooocooocrrrr ee 5 3 Geometry Algorithms 0 0 00 eee 5 4 PLOUFUNGCHIONS erie See A Pe ote 5 4 Utility Algorithms 0 0 cece eee 5 5 User Defined Algorithms 000 ccc eee eee 5 7 Demonstration ProgramS 00 cece eee ees 5 7 PDE Coefficients for Scalar Case 00005 5 20 PDE Coefficients for System Case 000000 5 21 Boundary Condition Dialog Box 20000 5 80 Mod Mfile co aa ett teases oes Wheels dia 5 81 Index Tutorial Introduction lt lt eee ee eee ee HAM eed A 1 2 What Does this Toolbox Do assassanin 1 2 Can I Use the PDE Toolbox 0 00 cee 1 2 What Problems Can Solve 0000 eee eee 1 3 In Which Areas Can the Toolbox Be Used 1 5 How Dol Definea PDE Problem 000 cee eee
32. 37 1 Tutorial Geometry Description matrix Decomposed Geometry matrix Geometry M file Boundary Condition matrix Mesh Ml refinemesh data assempde Boundary Coefficient Coefficient M file matrix M file Solution data pdeplot Constructive Solid Geometry Model A Constructive Solid Geometry CSG model is specified by a Geometry Description matrix a set formula and a Name Space matrix These data structures are described in the Reference chapter entry on dec sg At this level the problem geometry is defined by overlapping solid objects These can be created by drawing the CSG model in the GUI and then exporting the data using the Export Geometry Description Set Formula Labels option from the Draw menu 1 38 Using Command line Functions Decomposed Geometry A decomposed geometry is specified by either a Decomposed Geometry matrix or by a Geometry M file Here the geometry is described as a set of disjoint minimal regions bounded by boundary segments and border segments A Decomposed Geometry matrix can be created from a CSG model by using the function decsg It can also be exported from the GUI by selecting the Export Decomposed Geometry Boundary Cond s option from the Boundary menu A Geometry M file equivalent to a given Decomposed Geometry matrix can be created using the wgeo m function A decomposed geometry can be visualized with thepdegp ot function
33. Can use it What problems can solve etc What Does this Toolbox Do The Partial Differential Equation PDE Toolbox provides a powerful and flexible environment for the study and solution of partial differential equations in two space dimensions and time The equations are discretized by the Finite Element Method FEM The objectives of the PDE Toolbox are to provide you with tools that e Define a PDE problem i e define 2 D regions boundary conditions and PDE coefficients e Numerically solve the PDE problem i e generate unstructured meshes discretize the equations and produce an approximation to the solution e Visualize the results Can Use the PDE Toolbox The PDE Toolbox is designed for both beginners and advanced users The minimal requirement is that you can formulate a PDE problem on paper draw the domain write the boundary conditions and the PDE Start MATLAB At the MATLAB command line type pdetool This invokes the graphical user interface GUI which is a self contained graphical environment for PDE solving F or common applications you can use the specific physical terms rather than abstract coefficients Using pdet ool requires no knowledge of the mathematics behind the PDE the numerical schemes or MATLAB In Getting Started on page 1 9 we guide you through an example step by step Advanced applications are also possible by downloading the domain geometry boundary conditions and mes
34. Geometry Description Matrix Delete borders between minimal regions Decompose Constructive Solid Geometry into minimal regions Create initial triangular mesh J iggle internal points of a triangular mesh Interpolation between parametric representation and arc length Make regular mesh on a rectangular geometry Refine a triangular mesh Write boundary condition specification file Write geometry specification function Plot Functions Function Purpose pdecont pdegplot pdemesh pdepl ot pdesurf Shorthand command for contour plot Plot aPDE geometry Plot aPDE triangular mesh Generic PDE Toolbox plot function Shorthand command for surface plot Commands G rouped by Function Utility Algorithms Function Purpose dst Discrete sine transform idst Inverse discrete sine transform pdeadgsc Select triangles using a relative tolerance criterion pdeadworst pdecgrad pdeent pdegrad pdeintrp pdej mps pdeprtni pdesde pdesdp pdesdt pdesmech pdetrg Select triangles relative to the worst value The flux of a PDE solution Indices of triangles neighboring a given set of triangles The gradient of a PDE solution Interpolate from node data to triangle midpoint data Error estimates for adaption Interpolate from triangle midpoint data to node data Indices of edges in a set of subdomains Indices of points in a set of subdomains Indices of triangles in a set of subdom
35. M file ne pdegeom d pdegeom bs x y pdegeom bs s Werepresent 2 D regions by parameterized edge segments Both the regions and edge segments are assigned unique positive numbers as labels The edge segments cannot overlap The full 2 D problem description can contain several nonintersecting regions and they can have common border segments The boundary of a region can consist of several edge segments All edge segment junctions must coincide with edge segment endpoints We sometimes refer to an edge segment as a boundary segment or a border segment A boundary segment is located on the outer boundary of the union of the minimal regions and a border segment is located on the border between minimal regions There are two options for specifying the problem geometry e Create a Decomposed Geometry matrix with the function function decsg This is done automatically from pdet ool Using the Decomposed Geometry matrix restricts the edge segments to be straight lines circle or ellipse segments The Decomposed Geometry matrix can be used instead of the Geometry M file in the toolbox e Create a Geometry M file By creating your own Geometry M file you can create a geometry that follows any mathematical function exactly Below is an example of how to create a cardioid ne pdegeomis the number of edge segments d pdegeom bs isa matrix with one column for each edge segment specified in bs e Row 1 contains the start parameter value
36. PDE problem y can either be a Decomposed Geometry matrix or the name of a Geometry M file The formats of the Decomposed Geometry matrix and Geometry M file are described in the entries ondecsg and pdegeom respectively The outputs p e andt arethe mesh data In the Point matrixp the first and second rows contain x and y coordinates of the points in the mesh In the Edge matrix e the first and second rows contain indices of the starting and ending point the third and fourth rows contain the starting and ending parameter values the fifth row contains the edge segment number and the sixth and seventh row contain the left and right hand side subdomain numbers In the Triangle matrixt the first three rows contain indices to the corner points given in counter clockwise order and the fourth row contains the subdomain number inttmesh The following property name property value pairs are allowed Property Value Default Description Hma x numeric estimate Maximum edge size Hgrad numeric 1 3 Mesh growth rate Box on of f off Preserve bounding box Init on off off Edge triangulation Jiggle of f mean min mean Call j i ggl emesh Jigglelter numeric 10 Maximum iterations TheHmax property controls the size of the triangles on the mesh i ni t mesh creates a mesh where no triangle side exceeds H ma x Thedgr ad property determines te mesh growth rate away from a small part of the geometry The default value is 1 3 i
37. The snap to grid feature simplifies aligning the solid objects sj PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help aes AA AJS Generic Soatar X 00 Y 0 0 Grid Spacing lt lt Set forn Axes Limits Axes Equal Turn off Toolbar Help 0 Zoom Application e 0 Refresh 7 T i i i i i 15 1 0 5 0 05 1 15 Info Click and drag at perimeter to create ellipse Exit 1 The first step is to draw the geometry on which you want to solve the PDE The GUI provides four basic types of solid objects polygons rectangles circles and ellipses The objects are used to create a Constructive Solid Geometry model CSG model Each solid object is assigned a uniquelabel and by the use of set algebra the resulting geometry can be made up of a combination of unions intersections and set differences By default the resulting CSG model is the union of all solid objects 1 9 1 Tutorial To select a solid object either click on the button with an icon depicting the solid object that you want to use or select the object by using the Draw pull down menu In this case rectangle square objects are selected To draw a rectangle or a squarestarting at a corner press the rectangle button without a sign in the middle The button with the sign is used when you want to draw starting at the
38. The data structures of the Decomposed Geometry matrix and Geometry M file are described in Chapter 5 Reference entries on decsg andpdegeom respectively Boundary Conditions These are specified by either a Boundary Condition matrix or a Boundary M file Boundary conditions are given as functions on boundary segments A Boundary Condition matrix can be exported from the GUI by selecting the Export Decomposed Geometry Boundary Cond s option from the Boundary menu A Boundary M file equivalent to a given Boundary Condition matrix can be created using the wbound function The data structures of the Boundary Condition matrix and Boundary M file are described in Chapter 5 Reference entries onassemb and pdebound respectively Equation Coefficients The PDE is specified by either a Coefficient matrix or a Coefficient M file for each of the PDE coefficients c a f and d The coefficients are functions on the subdomains Coefficients can be exported from the GUI by selecting the Export PDE Coefficients option from the PDE menu The details on the equation coefficient data structures are given in the Chapter 5 Reference entry on assempde Mesh A triangular mesh is described by the mesh data which consists of a Point matrix an Edge matrix and a Triangle matrix In the mesh minimal regions are triangulated into subdomains and border segments and boundary segments are broken up into edges Mesh data is created from a decompos
39. This procedure yields a linear system KU F The components of U are the values of u at the nodes 1 20 Basics of The Finite Element M ethod For xinside a triangle u x is found by linear interpolation from the nodal values FEM techniques are also used to solve more general problems Below are some generalizations that you can access both through the graphical user interface and with command line functions e Time dependent problems are easy to implement in the FEM context The solution u x t of the equation EV eVu a f N can be approximated by u x t gt _ U x This yields a system of ordinary differential equations ODE Mou KU F which you integrate using ODE solvers Two time derivatives yield a second order ODE a TE dt derivatives the functions parabolic andhyperbolic M U KU F amp c The toolbox supports problems with one or two time Eigenvalue problems Solve V cVu au Adu for the unknowns u and A A is a complex number Using the FEM discretization you solve the algebraic eigenvalue problem KU 2 MU to find u and as approximations to u and A A robust eigenvalue solver is implemented in pdeei g If the coefficients c a f q or g are functions of u the PDE is called nonlinear and FEM yields a nonlinear system K U U F U You can use iterative methods for solving the nonlinear system The toolbox provides a nonlinear solver called pdenon in using a damped Gaus
40. a triangular mesh to a rectangular grid pdegrad and pdecgrad compute gradients of the solution pdepl ot has a large number of options for plotting the solution pdecont and pdesurf are convenient shorthands for pdep ot Hints and Suggestions for Using Command Line Functions Several examples of command line function usage are given in the Examples chapter Use the export facilities of the GUI as much as you can They provide data structures with the correct syntax and these are good starting points that you can modify to suit your needs A good way to produce a Geometry M file describing a geometry outside of the possibilities provided by the GUI is to draw a similar geometry using the GUI export the Decomposed Geometry matrix and write a Geometry M file with wgeom The special segments can then be edited by hand An example of a 1 40 Using Command line Functions hand tailored Geometry M fileis car dg also presented in the Chapter 5 Reference entry on pdegeom Working with the system matrices and vectors produced by as s e ma and assemb can sometimes be valuable When solving the same equation for different loads or boundary conditions it pays to assemble the stiffness matrix only once Point loads on a particular node can be implemented by adding the load tothe corresponding row in the right hand side vector A non local constraint can be incorporated intotheH and R matrices An example of a hand written Coefficient
41. an animation check the Animation check box in the Plot selection dialog box Also select the colormap hot Press the Plot button to start a recording of the solution plots in a separate figure window The recorded animation is then played five times Note that the temperature in the block rises very quickly To improve the animation and focus on the first second try to change the list of times to the MATLAB expression ogspace 2 0 5 20 2 18 Examples of Parabolic Problems Also try to change the heat capacity coefficient d and the heat flow at the right most boundary to see how they affect the heat distribution PDE Toolbox Untitled Time 2 5 Color u TTT x i i iii i i 15 1 0 5 0 09 05 Os 1 15 Info Select a different plot or change mode to alter PDE mesh or boundary conditions Exit Using Command Line Functions First you must create geometry and boundary condition M files The M files used here were created using pdet ool The geometry of the metal block is described in crackg m and the boundary conditions can be found incrackb m To create an initial mesh call i ni t mesh p e t initmesh crackg The heat equation can now be solved using the PDE Toolbox function parabolic The generic parabolic PDE that parabolic solves is du de 7 V cVu au f with initial value u0 u t and the times at which to compute a sol
42. buttons and they represent from left to right e Draw a rectangle square starting at a corner e Draw a rectangle square starting at the center e Draw an ellipse circle starting at the perimeter e Draw an ellipse circle starting at the center e Draw a polygon Click and drag to create polygon sides You can close the polygon by pressing the right mouse button Clicking at the starting vertex also closes the polygon The Draw mode buttons can only be activated one at the time and they all work the same way single clicking on a button allows you to draw one solid object of the selected type Double clicking on a button makes it stick and you can then continue to draw solid objects of the selected type until you single click on the button to release it Using the right mouse button or Control click the drawing is constrained to a square or a circle 1 25 1 Tutorial The second group of six buttons includes the following buttons e Q Enters the Boundary mode PDE Opens the PDE Specification dialog box e A Initializes the triangular mesh e The refine button the four triangle icon Refines the triangular mesh e Solves the PDE e The 3 D solution plot icon Opens the Plot Selection dialog box The right most button with the magnifying glass toggles the zoom function on off The GUI Modes The PDE solving process can be divided into several steps 1 Define the geometry 2 D domain 2 Define t
43. example of an elliptic problem let s use the simplest elliptic PDE of all Poisson s equation The problem formulation is AU 1inQ U 00n dQ where Q is the unit disk In this case the exact solution is 1 x y U x y 7 so the error of the numeric solution can be evaluated for different meshes Using the Graphical User Interface With thepdetoo GUI started perform the following steps using the Generic Scalar mode 1 Using some of the Option menu features add a grid and turn on the snap to grid feature Draw a circle by clicking on the button with the ellipse icon with the sign and then click and drag from the origin using the right mouse button to a point at the circle s perimeter If the circle that you createis not a perfect unit circle double click the circle This opens a dialog box where you can specify the exact center location and radius of the circle 2 Enter the Boundary mode by clicking on the button with the Q icon The boundaries of the decomposed geometry are plotted and the outer boundaries are assigned a default boundary condition Dirichlet boundary condition u 0 on the boundary In this case this is what we want If the boundary condition is different double click on the boundary to open a dialog box through which you can enter and display the boundary condition Examples of Elliptic Problems Todefinethe partial differential equation press the PDE button This opens a dialog
44. file q g h r pdebound p e u ti me computes the values of q g h and r on the a set of edges e The matricesp ande are mesh data e needs only to be a subset of the edges in the mesh Details on the mesh data representation can be found in the entry oninitmesh The input arguments u andti me are used for the nonlinear solver and time stepping algorithms respectively u and ti me are empty matrices if the corresponding parameter is not passed toassemb Ifti me isNaN and any of the function q g h and r depends ont i me pdebound must return a matrix of correct size containing NaNsin all positions in the corresponding output argument The solution u is represented by the solution vector u Details on the representation can be found in the entry onassempde pdebound Examples See Also q and g must contain the value of q and g on the midpoint of each boundary Thus we havesize q N 2 ne whereN isthe dimension of the system and ne the number of edges ine andsize g N ne For the Dirichlet case the corresponding values must be zeros h andr must contain the values of h and r at the first point on each edge followed by the value at the second point on each edge Thus we have size h N 2 2 ne wheren is the dimension of the system and ne the number of edgesine andsize r N 2 ne When M lt N h and r must be padded with N M rows of zeros The elements of the matrices q and h are stored in column wise ordering
45. fixed point iteration is also possible as follows Essentially for a given U compute the FEM matrices K and F and set yard K IF This is equivalent to approximating the J acobian with the stiffness matrix Indeed since p U KU F putting y K yields n 1 U DUO ERA In many cases the convergence rate is slow but the cost of each iteration is cheap 4 23 4 The Finite Element M ethod The nonlinear solver implemented in the PDE Toolbox also provides for a compromise between the two extremes To compute the derivative of the mapping U gt KU proceed as follows The a term has been omitted for clarity but appears again in the final result below o0 KU 1 f U VO VO dxU a c U VO VO oy E Vo Vo dxU The first integral term is nothing more than K ij The second term is lumped i e replaced by a diagonal matrix that contains the row sums Since 250 1 the second term is approximated by DN 5 V9 Vo dxU which is the itcomponent of KU where K is the stiffness matrix associated with the coefficient ge rather than c The same reasoning can be applied to the derivative of the mapping U gt MU Finally note that the derivative of the mapping U gt F is exactly j 2 Lo 9jdx of which is the mass matrix associated with the coefficient Sy Thus the acobian of the residual p U is approximated by J K MCP y diags KO mM 0 where the differentiation is
46. for modeling and solving your PDE problems using thepdetool GUI Therearealsoa number of shortcuts that you can use in certain situations NOTE There are platform dependent keyboard accelerators available for many of the most common pdet oo GUI activities Learning to use the accelerator keys may improve the efficiency of your pdetool sessions The basic flow of actions is indicated by the way the graphical push buttons and the menus are ordered from left to right You work your way from left to right in the process of modeling defining and solving your PDE problem using the pdetool GUI e When you start pdetool isin a Draw mode where you can use the four basic solid objects to draw your Constructive Solid Geometry CSG model Y ou can also edit the set formula The solid objects are selected using the five left most push buttons or from the Draw menu e Totheright of the Draw mode buttons you find push buttons through which you can access all the functions that you need to define and solve the PDE problem define boundary conditions design the triangular mesh solve the PDE and plot the solution The following sequence of actions covers all the steps of a normal pdet ool session 1 Usepdetool asa drawing tool to make a drawing of the 2 D geometry on which you want to solve your PDE Make use of the four basic solid objects and the grid and the snap to grid feature The GUI starts in the Draw mode and you can select
47. for the L Shaped Membrane The problem of finding the eigenvalues and the corresponding eigenfunctions of an L shaped membrane is of interest to all MATLAB users since the plot of the first eigenfunction is the logo of The MathWorks In fact you can compare the PDE Toolbox computed eigenvalues and eigenfunctions to the ones produced by the MATLAB function membrane The problem is to compute all eigenmodes with eigenvalues lt 100 for the agenmode PDE problem Au Au on the geometry of the L shaped membrane u 0 on the boundary Dirichlet condition Using the Graphical User Interface With thepdetoo l GUI active check that the current mode is Generic Scalar Then draw the L shape as a polygon with corners in 0 0 1 0 1 1 1 1 1 1 and 0 1 There is no need to define any boundary conditions for this problem since the default condition u 0 on the boundary is the correct one Therefore you can continue to the next step to initialize the mesh Refine the initial mesh twice Defining the eigenvalue PDE problem is also easy Open the PDE Specification dialog box and select Eigenmodes The default values for the PDE coefficients c 1 a 0 d 1 all match the problem description so you can exit the PDE Specification dialog box by pressing the OK button Open the Solve Parameters dialog box by selecting Parameters fromthe Solve menu The dialog box contains an edit box for entering the eigenvalue search ra
48. function of the solution u you must use the nonlinear solver If the boundary condition is a function of the timet you must choose a parabolic or hyperbolic PDE Examples 100 80 s nx andcos x 2 In the nongeneric application modes the Description column contains descriptions of the physical interpretation of the boundary condition parameters 3 17 3 The Graphical User Interface 3 18 PDE Menu PDE Mesh Solve Plot PDE Mode Show Subdomain Labels PDE Specification Export PDE Coefficients PDE Mode Show Subdomain Labels PDE Specification Export PDE Coefficients Enter the partial differential equation mode Toggle the labeling of the subdomains on off The subdomains are labeled using the subdomain numbering in the decomposed geometry matrix Open dialog box for entering PDE coefficients and types Export current PDE coefficients tothe main workspace The resulting workspace variables are strings PDE Toolbox M enus PDE Specification 5 PDE Specification Equation div c grad u a u f Type of PDE Coefficient Value 4 Elliptic c 1 0 wv Parabolic a 0 0 w Hyperbolic f 10 v Eigenmodes i OK Cancel PDE Specification opens a dialog box where you enter the type of partial differential equation and the applicable parameters The dimension of the parameters is dependent on the dimension of the PDE The description below a
49. icon or select Parameters from the Plot menu to access the dialog box for selection of the different plot options Several plot styles are available and the solution can be plotted in the GUI or in a separate figure as a 3 D plot Now select a plot where the color and the height both represent u Choose interpolated shading and use the continuous interpolated height option The default colormap is thecoo colormap a pop up menu lets you select from a number of different colormaps Finally press the Plot button to plot the solution press the Done button to save the plot setup as the current default The solution is plotted as a 3 D plotina separate figure window G etting Started y Plot Selection Plot type Property User entry Plot style E Color r AX interpolated shad Contour _ Arrows grad u al proportional 3 Deformed mesh grad u Height 3 D plot u continuous 4 at Plot in x y grid Contour plot levels 20 E Plot solution automatically Show mesh Colormap cool i Plot Done Cancel al The following solution plot is the result You can use the mouse to rotate the plot in 3 D By clicking and dragging the axes the angle from which the solution is viewed can be changed y Figure No 3 Color u Height u This concludes the first example of solving a PDE by using thepdetoo GUI Many more examples in Chapter 2 Examples focus on sol
50. in the MATLAB matrices q and h For the boundary conditions 1 1 u 2 2 a 3 Jew 20 4 the values below should be stored in q g h andr 1 q Ze 2 0 g 33 4 1 1 h a 0 1 1 0 0 r oe D3 2 0 0 initmesh pdegeom pdesdt pdeent pdecgrad Purpose Synopsis Description See Also 5 50 The flux of a PDE solution cgxu cgyu pdecgrad p t c u cgxu cgyu pdecgrad p t c u ti me cgxu cgyu pdecgrad p t c u time sdl cgxu cgyu pdecgrad p t c u returns the flux c O Vu evaluated at the center of each triangle Row i of cgxu contains N u du PE ill Fx y j l Row i of cgyu contains N du du L Ej ij29y j 1 There is one column for each triangleint in both cgxu andcgyu The geometry of the PDE problem is given by the mesh data p andt Details on the mesh data representation can be found in the entry oni ni t mesh The coefficient of the PDE problem can be given in a variety of ways A complete listing of all options can be found in the entry on assempde The format for the solution vector u is described inassempde The scalar optional argument ti me is used for parabolic and hyperbolic problems ifc depends on t the time The optional argument sd restricts the computation to the subdomains in the list sdi assempde pdecirc Purpose Synopsis Description Examples See Also Draw circle pdecirc xc yc radius pd
51. length N The elements Cijkl Aij dij MATLAB matricesc a d andf Thecase of identity diagonal and symmetric matrices are handled as special cases For the tensor Gj this applies both to the indices i andj and tothe indices k and I and f of c a d and f are stored row wise in the The PDE Toolbox does not check the ellipticity of the problem and it is quite possible to define a system that is not elliptic in the mathematical sense The procedure described above for the scalar case is applied to each component of the system yielding a symmetric positive definite system of equations whenever the differential system possesses these characteristics 4 10 The Elliptic System The boundary conditions now in general are mixed e for each point on the boundary a combination of Dirichlet and generalized Neumann conditions hu r c Vu qu g hk u By the notation c Vu we mean the N by 1 matrix with i 1 component N lo 0 lo lo y cos Q c 1 15 cos O c i i sin O c 213 sin A c 2 23 j 1 where the outward normal vector of the boundary is cos Q sin a There are M Dirichlet conditions and the h matrix is M by N M gt 0 The generalized Neumann condition contains a source h u where the Lagrange multipliers u are computed such that the Dirichlet conditions become satisfied In a structural mechanics problem this term is exactly the reaction force neces
52. of 0 2 in the center and by outer limits with side length of 0 5 At the inner boundary the electrostatic potential is LOOOV At the outer boundary the electrostatic potential is OV Thereis no charge in the domain This leads to the problem of solving the Laplace equation AV 0 with the Dirichlet boundary conditions V 1000 at the inner boundary and V O at the outer boundary Using the Graphical User Interface After selecting the application mode Electrostatics the 2 D area is most easily drawn by first drawing a square with sides of length 0 2 use the Snap option and adjust the grid spacing if necessary Then draw another square with sides of length 0 5 using the same center position The 2 D domain is then simply SQ2 SQ1 if the first square is named SQ1 and the second square is named SQ2 Enter the expression into the Set formula edit box and proceed to define the boundary conditions Use Shift click to select all the inner boundaries Then double click on an inner boundary and enter 1000 as the Dirichlet boundary condition for the inner boundaries Next open the PDE Specification dialog box and enter 0 into the space charge density r ho edit field The coefficient of dielectricity can be left at 1 since it does not affect the result as long as it is constant Initialize the mesh and press the button to solve the equation Using the adaptive mode you can improve the accuracy of the solution by refining the mesh close t
53. of the eigenvalues of this Hessenberg matrix Hj j eventually give good approximations to the eigenvalues of the original pencil K M when the basis grows in dimension j and less and less of the eigenvector is hidden in the residual matrix Ej The basis V is built one column vj at a time The first vector v is chosen at random as n normally distributed random numbers In stepj thefirstj vectors are already computed and form the n xj matrix Vj The next vector vj 1 is computed by first letting A operate on the newest vector vj and then making the result orthogonal to all the previous vectors This is formulated as hj 1jVj 1 Av Vjhj where the column vector hj consists of the Gram Schmidt coefficients and hj 41 is the normalization factor that gives vj 1 unit length Put the corresponding relations from previous steps in front of this and get T h AV VH y e V V iY j 1 jfj j 1 where Hj j is aj xj Hessenberg matrix with the vectors hj as columns The second term on the right hand side has nonzeros only in the last column the earlier normalization factors show up in the subdiagonal of Hj j The eigensolution of the small Hessenberg matrix Hj j gives approximations to some of the eigenvalues and eigenvectors of the large matrix operator A in the following way Compute eigenvalues 6 and eigenvectors sj of Hj j 4 18 The Eigenvalue Equation HS S9 i 1 eta ae Then y Vjsjis an approximate eigenvector of
54. of the nodal values If the mesh is uniformly refined Vy approximates the set of smooth functions on Q A suitable basis for Ya is theset of tent or hat functions o These are linear on each triangle and take the value 0 at all nodes x except for xj Requesting 0 x 1 yields the very pleasant property The Elliptic Equation N u x y Uo x U j 1 i e by solving the FEM system we obtain the nodal values of the approximate solution Finally note that the basis function q vanishes on all the triangles that do not contain the node x The immediate consequence is that the integrals appearing in K ij Mij Qij Fj and G only need to be computed on the triangles that contain the node x Secondly it means that K j and M j arezero unless x and x are vertices of the same triangle and thus K and M are very sparse matrices Their sparse structure depends on the ordering of the indices of the mesh points The integrals in the FEM matrices are computed by adding the contributions from each triangle to the corresponding entries i e only if the corresponding mesh point is a vertex of the triangle This process is commonly called assembling hence the name of the function assempde The assembling routines scan the triangles of the mesh F or each triangle they compute the so called local matrices and add their components to the correct positions in the sparse matrices or vectors The integrals are computed using
55. or by selecting PDE Specification fromthePDE menu This brings up a dialog box In Mesh mode you can control the automated mesh generation and plot the mesh An initial mesh can be generated by clicking on the A button or by selecting Initialize Mesh from the Mesh menu The initial mesh can be repeatedly refined by clicking on the refine button or by selecting Refine Mesh from the Mesh menu In Solve mode you can specify solve parameters and solve the PDE For parabolic and hyperbolic PDE problems you can also specify the initial conditions and the times at which the output should be generated F or eigenvalue problems the search range can be specified Also the adaptive and nonlinear solvers for elliptic PDEs can be invoked The PDE problem is solved by clicking on the button or by selecting Solve PDE from the Solve menu By default the solution is plotted in thepdetool axes 5 79 pdetool 5 80 In Plot mode you can select a wide variety of visualization methods such as surface mesh contour and quiver vector field plots For surface plots you can choose between interpolated and flat rendering schemes The mesh can be hidden in all plot types For parabolic and hyperbolic equations you can animate the solution as it changes with time You can show the solution both in 2 D and 3 D 2 D plots are shown insidepdet oo 3 D plots are plotted in separate figure windows Different types of plots can be selected by clicking on
56. orthogonalized against all the 22 converged vectors the smallest eigenvaluestabilized at a value outside of theinterval 0 100 sothe algorithm signaled convergence Of the 22 converged eigenvalues 19 were inside the search interval L Shaped Membrane with Rounded Corner An extension of this problem is to compute the eigenvalues for an L shaped membrane where the inner corner at the knee is rounded The roundness is created by adding a circle so that the circle s arc is a part of the L shaped membrane s boundary By varying the circle s radius the degree of roundness can be controlled The M filel shapec is an extension of an ordinary model file created using pdet ool It contains the lines pdepoly 1 1 1 0 0 1 Pel L Ty L 0 OP Pas pdecirc a a a Cl pderect a 0 a 0 SQ1 The extra circle and rectangle that are added using pdecirc andpderect to create the rounded corner are affected by the added input argument a through a couple of extra lines of MATLAB code This is possible since the PDE Toolbox is a part of the open MATLAB environment With shapec you can create L shaped rounded geometries with different degrees of roundness If you usel shapec without an input argument a default radius of 0 5 is used Otherwise use shapec a wherea is the radius of the circle Experimenting using different values for the radius a shows you that the eigenvalues and the frequencies of the corresponding eig
57. pdesmech 5 75 pdesurf 5 77 pdetool 5 78 pdetrg 5 83 pdetriq 5 84 pencil 5 92 plane strain 2 41 plane stress 2 36 Plot menu 3 30 plot mode 1 27 Plot parameters 3 30 poiasma 5 85 poicalc 5 86 poiindex 5 87 poimesh 5 88 Point matrix 1 39 5 38 poisolv 5 89 Poisson s equation 4 27 Poisson s equation 1 9 2 2 2 55 polygon solid 5 31 Print command 3 6 R rectangle solid 5 32 refinemesh 5 90 Robin boundary condition 4 3 Rotate command 3 14 S Save As command 3 5 scattering problem 2 6 Schur vector 4 19 set formula 1 5 1 11 1 28 1 38 5 32 skin effect 2 52 solid object 1 9 1 27 5 31 5 78 solution vector 1 40 5 19 Solve menu 3 25 solve mode 1 27 solve parameters 3 25 1 3 Index 1 4 sptarn 5 92 stiff springs 4 11 5 19 stiffness matrix 4 6 4 16 structural mechanics 2 36 2 41 T test function 1 19 4 4 toolbar 1 25 3 38 tri2grid 5 95 triangle data 5 62 5 70 Triangle matrix 1 39 5 38 triangle quality 3 22 5 84 V von Mises effective stress 2 39 2 42 5 76 W wave equation 2 23 wbound 5 96 weak form 4 4 wgeom 5 97 Window menu 3 37
58. polygon is provided You can combine these objects using set formulas e n Boundary mode you specify the boundary conditions You can have different types of boundary conditions on different boundary segments e In PDE mode you interactively specify the type of PDE and the coefficients c a f and d You can specify the coefficients for each subdomain independently This may ease the specification of e g various material properties in a PDE model 1 Tutorial How Can I Solve a PDE Problem Most problems can be solved from the graphical user interface There are two major modes that help you solve a problem e n Mesh mode you generate and plot meshes You can control the parameters of the automated mesh generator e n Solve mode you can invoke and control the nonlinear and adaptive sol vers for elliptic problems For parabolic and hyperbolic problems you can specify the initial values and the times for which the output should be generated For the eigenvalue solver you can specify the interval in which to search for eigenvalues After solving a problem you can return to the Mesh mode to further refine your mesh and then solve again Y ou can also employ the adaptive mesh refiner and solver This option tries to find a mesh that fits the solution Can I Use the Toolbox for Nonstandard Problems For advanced nonstandard applications you can transfer the description of domains boundary conditions etc to your MATLAB
59. r v x y e In this problem the incident wave is traveling in the x direction so the boundary condition is simply r e Enter this boundary condition in the Boundary Condition dialog box as a Dirichlet condition h 1 1 exp i 60 x Thereal part of this is a sinusoid For sufficient accuracy about 10 finite elements per wavelength are needed The outer boundary should be located a few object diameters from the object itself An initial mesh generation and two successive mesh refinements give approximately the desired resolution Although originally a wave equation the transformation into a Helmholtz s equation makes it in the PDE Toolbox context but not strictly mathematically an elliptic equation The elliptic PDE coefficients for this problem arec 1 a k 3600 and f 0 Open the PDE Specification dialog box and enter these values The problem can now be solved and the solution is complex For a complex solution the real part is plotted and a warning message is issued The propagation of the reflected waves is computed as Re r x y e which is the reflex of 2 8 Examples of Elliptic Problems i ka ot e Rel Toseethe whole field plot ika iot Re r x y e je The reflected waves and the shadow behind the object are clearly visible when you plot the reflected wave Tomakean animation of the reflected wave the solution and the mesh data must first be exported to
60. should be generated For the eigenvalue solver you can specify the interval in which to search for eigenvalues In Plot mode there is wide range of visualization possibilities You can visualize both in thepdet 00 GUI and in a separate figure window You can visualize three different solution properties at the same time using color height and vector field plots There are surface mesh contour and arrow quiver plots available For parabolic and hyperbolic equations you can animate the solution as it changes with time The CSG Model and the Set Formula ThePDE Toolbox uses the Constructive Solid Geometry CSG model paradigm for the modeling You can draw solid objects that can overlap There are four types of solid objects e Circle object represents the set of points inside and on a circle e Polygon object represents the set of points inside and on a polygon given by a set of line segments e Rectangle object represents the set of points inside and on a rectangle e Ellipse object represents the set of points inside and on an ellipse The ellipse can be rotated Each solid object is automatically given a unique name by the GUI The default names are C1 C2 C3 etc for circles P1 P2 P3 etc for polygons R1 R2 R3 etc for rectangles E1 E2 E3 etc for ellipses Squares although just a special case of rectangles are named SQ1 SQ2 SQ3 etc The nameis displayed on the solid object itself You c
61. so repeated Undo Mesh Change eventually bring you back to the initial mesh Display a plot of the triangular mesh where the individual triangles are colored according to their quality The quality measure is a number between 0 and 1 where triangles with a quality measure greater than 0 6 are acceptable F or details on the triangle quality measure see pdetriq in the Reference chapter PDE Toolbox M enus Show Node Labels Show Triangle Labels Parameters Export Mesh Parameters esh Parameters gt MeshParameters Initmesh parameters Maximum edge size E Jiggle mesh Jigglemesh parameters Jiggle mode optimize mean Number of jiggle iterations Refinement method regular ok Cancel el Toggle the mesh node labels on off The node labels are the column numbers in the Point matrix p Toggle the mesh triangle labels on off The triangle labels are the column numbers in the triangle matrix t Open dialog box for modification of mesh generation parameters Export Point matrix p Edge matrix e and Triangle matrix t to the main workspace Parameters opens a dialog box containing mesh generation parameters The parameters used by the mesh initialization algorithm ni t mesh are e Maximum edge size Largest triangle edge length approximately This parameter is optional and must be a real positive number e Mesh growth r
62. t c a f u alfa beta m errf pdej mps p t c a f u alfa beta m calculates the error indication function used for adaption The columns of er rf correspond to triangles and the rows correspond to the different equations in the PDE system p and t are mesh data Seethe entry oni ni t mesh for details c a and f are PDE coefficients See the entry on assempde for details c a andf must be expanded so that columns correspond to triangles u is the solution vector See the entry on assempde for details The formula for computing the error indicator E K for each triangle K is 1 2 EK of alet B Y helm Yup te OK where n is the unit normal of edge t and the braced term is the jump in flux across the element edge The norm is an Lz norm computed over the element K Theerror indicator is stored in er rf as column vectors one for each triangle int Moreinformation can be found in the section Adaptive M esh Refinement on page 4 26 adapt mesh pdeadgsc pdeadworst 5 63 pdemdlcv Purpose Synopsis Description Example 5 64 Convert PDE Toolbox 1 0 Model M files to PDE Toolbox 1 0 2 format pdemdIcv infile outfile pdemdil cv infile outfile converts the PDE Toolbox 1 0 Model fileinfi le toa PDE Toolbox 1 0 2 compatible Model M file The resulting M file is saved asoutfile Ifthe m extension is missing in out file itis added automatically If you want to use Model M files generated using PDE Toolbox 1
63. the GUI The edit box for the set formula contains the active set formula In the main axes you draw the 2 D geometry display the mesh plot the solution etc At the bottom of the GUI an information line provides information about the current activity It can also display help information about the toolbar buttons The Menus There are 11 different pull down menus in the GUI For a more detailed description of the menus and the dialog boxes see Chapter 3 The Graphical User Interface File menu From the File menu you can Open and Save model M files that contain a command sequence that reproduces your modeling session Y ou can also print the current graphics and exit the GUI Edit menu From the E dit menu you can cut clear copy and paste the solid objects There is also a Select All option Options menu The Options menu contains options such as toggling the axis grid a snap to grid feature and zoom Y ou can also adjust the axis limits and the grid spacing select the application mode and refresh the GUI Draw menu From the Draw menu you can select the basic solid objects such as circles and polygons You can then draw objects of the selected type using the mouse From the Draw menu you can also rotate the solid objects and export the geometry to the MATLAB main workspace Boundary menu From the Boundary menu you access a dialog box where you define the boundary conditions Additionally you can label edges and subdo
64. thei ni t mesh function The three components of the algorithm are the error indicator function which computes an estimate of the element error contribution the mesh refiner which selects and subdivides elements and the termination criteria The Error Indicator Function The adaption is a feedback process As such it is easily applied to a larger range of problems than those for which its design was tailored You want estimates selection criteria etc to be optimal in the sense of giving the most accurate solution at fixed cost or lowest computational effort for a given accuracy Such results have been proved only for model problems but generally the equidistribution heuristic has been found near optimal Element sizes should be chosen such that each element contributes the same to the error The theory of adaptive schemes makes use of a priori bounds for Adaptive M esh Refinement solutions in terms of the source function f For nonelliptic problems such a bound may not exist while the refinement scheme is still well defined and has been found to work well The error indicator function used in the toolbox is an element wise estimate of the contribution based on the work of C J ohnson et al 5 6 For Poisson s equation Au f on Q the following error estimate for the FEM solution up holds in theL norm Viu u lt 0 hfi BD u where h h x is the local mesh size and D v Y m109 an 17 TEE
65. usingn levels pdecont p t u v plots using the levels specified by v This command is just shorthand for the call pdeplot p t xydata u xystyle off contour on levels n colorbar off If you want to have more control over your contour plot usepdep ot instead of pdecont Plot the contours of the solution to the equation Au 1 over the geometry defined by the L shaped membrane Use Dirichlet boundary conditions u 0 on Q p e t initmesh Ishapeg p e t refinemesh Ishapeg p e t usassempde shapeb p e t 1 0 1 pdecont p t u pdeplot pdemesh pdesurf pdeeig Purpose Synopsis Description Solve eigenvalue PDE problem v l pdeeig b p e t c a d r v pdeeig K B M r v pdeeig b p e t c a d r produces the solution tothe FEM formulation of the scalar PDE eigenvalue problem A cVu au Adu on Q or the system PDE eigenvalue problem V c 8 Vu au Adu on Q on a geometry described by p e andt and with boundary conditions given byb r isatwoelement vector indicating an interval on the real axis The left hand side can be nf The algorithm returns all eigenvalues in this interval in v iS an genvector matrix For the scalar case each column in v is an eigenvector of solution values at the corresponding node points fromp Fora system of dimension N with np node points the first np rows of v describe the first component of
66. ut0 tlist Kk F B ud M rtol ul hyperbolic u0 ut0 tlist Kk F B ud M rtol atol ul hyperbolic u0 ut0 tlist b p e t c a f d produces thesolution tothe FEM formulation of the scalar PDE problem d V cVu tau f on Q or the system PDE problem qo V C9 Vu tau f ong ar on a mesh described by p e and t with boundary conditions given by b and with initial valueuo and initial derivative ut 0 In the scalar case each row in the solution matrix u1 is the solution at the coordinates given by the corresponding column in p Each column in u1 is the solution at the time given by the corresponding itemint1Iist For a system of dimension N with np node points the first np rows of u1 describe the first component of u the following np rows of u1 describe the second component of u and so on Thus the components of u are placed in blocks u as N blocks of node point rows b describes the boundary conditions of the PDE problem b can be either a Boundary Condition matrix or the name of a Boundary M file The boundary conditions can depend ont the time The formats of the Boundary Condition matrix and Boundary M file are described in the entries on assemb and pdebound respectively The geometry of the PDE problem is given by the mesh datap e andt Details on the mesh data representation can be found in the entry oni nit mesh The coefficientsc a d andf of the PDE problem can be given in a variety of ways The coefficients ca
67. v the following np rows of v describe the second component of v and so on Thus the components of v are placed in blocks y as N blocks of node point rows b describes the boundary conditions of the PDE problem p can be either a Boundary Condition matrix or the name of a Boundary M file The formats of the Boundary Condition matrix and Boundary M file are described in the entries on assemb and pdebound respectively Note that the eigenvalue PDE problem is a homogeneous problem i e only boundary conditions where g 0 and r 0 can be used The non homogeneous part is removed automatically The geometry of the PDE problem is given by the mesh datap e andt Details on the mesh data representation can be found in the entry oni nit mesh The coefficients c a d of the PDE problem can be given in a wide variety of ways In the context of pdeei g the coefficients cannot depend on u nor t the time A complete listing of all options can be found in the entry on assempde v pdeeig K B M r produces the solution to the generalized sparse matrix eigenvalue problem K u AB MBu u Bu with Real A in the interval in r 5 53 pdeeig Examples Cautionary See Also 5 54 Compute the eigenvalues less than 100 and corresponding eigenmodes for Vu Au on the geometry of the L shaped membrane Then display the first and sixteenth eigenmodes p e t initmesh Ishapeg p e t refinemesh Ilshapeg p e t
68. vector of length si ze p 2 and triangle data is represented by a row vector of length size t 2 Ifnoxydata zdata orf owdata is Supplied pdepl ot plots the mesh specified by p e andt The option mesh displays or hides default the triangle mesh in the plot The option xy grid first converts the data to x y data usingtri2grid and then uses a standard MATLAB plotting algorithm The property gri dpar am passes thetri2grid data topdepl ot This speeds up animation see pdedemo5 and pdedemo6 The property col or map renders the plot using any MATLAB colormap or color matrix col orbar adds a colorbar tothe plot ti tle inserts a title into the plot evel s only applies to contour plots Given a scalar integer value it plots that number of equally spaced contour levels given a vector of level values it plots that contour lines on the levels in the vector h pdepl ot p t u additionally returns handles to the drawn axes objects The command sequence below plots the solution to Poisson s equation on the L shaped membrane in 3 D p e t initmesh Ishapeg u assempde Ishapeb p e t 1 0 1 pdeplot p e t xydata u zdata u mesh off pdecont pdegplot pdemesh pdesurf pdepoly Purpose Synopsis Description Examples See Also Draw polygon pdepoly x y pdepoly x y label pdepoly x y draws a polygon with corner coordinates defined by x andy If thepdetool GUI is not active it is automatica
69. with N rows one row for each component The case of identity diagonal and symmetric matrices are handled as special cases For the tensor Gijk this applies both to the indices i and j and to the indices k and l Thenumber of rowsinf determines the dimension N of the system Row i inf represents the component fi in f The number of rows ng ina is related to the components ajj of a according to the following table For the symmetric case assume that j gt i All elements aj that cannot be formed are assumed to be Zero Na Symmetric aij Row in a 1 No aij 1 N No aij i N N 1 2 Yes aij jG D 2 i N No aij N j 1 i 5 21 assempde 5 22 An example of how a is stored ina can be found below The coding of c in c is determined by the dimension N and the number of rows n inc The number of rows n in c is matched to the function of N in the first column in the table below sequentially from the first line to the last line The first match determines thetype of coding of c This actually means that for some small values 2 lt N lt 4 the coding is only determined by the order in which the tests are performed For the symmetric case assume that j gt i and gt k All elements cj that can not be formed are assumed to be zero Ne Symmetric Giki Row in c 1 No Gikk 1 2 No Gikk k 3 Yes Ciki I k 1 4 No Giikl 2l k 2 N No Gikk i 2N No Gikk 2i k 2 3N Yes Giikl 3i 1 k 4 4N No Gik 4i 21 k 6 2
70. workspace F rom there you use the functions of the PDE Toolbox for managing data on unstructured meshes You have full access to the mesh generators FEM discretizations of the PDE and boundary conditions interpolation functions etc You can design your own solvers or use FEM tosolvesubproblems of more complex algorithms See also the section Using Command Line F unctions How Can I Visualize My Results From the graphical user interface you can use Plot mode where you havea wide range of visualization possibilities You can visualize both inside the pdetool GUI and in separate figures You can plot three different solution properties at the same time using color height and vector field plots Surface mesh contour and arrow quiver plots are available For surface plots you can choose between interpolated and flat rendering schemes The mesh may be hidden or exposed in all plot types F or parabolic and hyperbolic equations you can even produce an animated movie of the solution s time dependence All visualization functions are also accessible from the command line 1 6 Intro duction Are There Any Applications Already Implemented The PDE Toolbox is easy to use in the most common areas due to the application interfaces Eight application interfaces are available in addition to the generic scalar and system vector valued u cases e Structural Mechanics Plane Stress e Structural Mechanics Plane Strain e Electr
71. y direction Look at the difference between the first and the second eigenvalue 1 2 1 1 ans 2 4740 pi pis 4 ans 2 4674 Likewise the fifth eigenmode is made up of the first eigenmode in the x direction and the third eigenmode in the y direction As expected 5 1 1 is approximately equal to 1 You can explore higher modes by increasing the search range to include eigenvalues greater than 10 2 34 Application Modes Application Modes In this section we describe the application modes of the pdet ool GUI Examples are given for a variety of applications to different engineering problems The Application Modes and the GUI The PDE Toolbox can be applied to a great number of problems in engineering and science To help you in using thepdet oo GUI for some important and common applications eight different application modes are available in addition to the generic scalar and system modes The available application modes are e Generic scalar the default mode e Generic system e Structural Mechanics Plane Stress Structural Mechanics Plane Strain e Electrostatics M agnetostatics e AC Power Electromagnetics e Conductive Media DC Heat Transfer Diffusion NOTE From the GUI the system PDEs are restricted to problems with vector valued u of dimension two Using command line functions there is no formal restriction on the dimension of u The application mode can be selected directly fro
72. 0 they must first be converted using pdemdl cv pdedmdl cv model 42 m model 5 m converts the PDE Toolbox 1 0 Model M file model 42 m and saves the converted model in model 5 m pdemesh Purpose Synopsis Description Examples See Also Plot a PDE triangular mesh pdemesh p e t pdemesh p e t u h pdemesh p e t h pdemesh p e t u pdemesh p e t plots the mesh specified by the mesh data p e andt h pdemesh p e t additionally returns handles to the plotted axes objects pdemesh p e t u plots PDE node or triangle data u using a mesh plot Ifu is a column vector node data is assumed If u is a row vector triangle data is assumed This command plots substantially faster than thepdes ur f command The geometry of the PDE problem is given by the mesh datap e andt Details on the mesh data representation can be found in the entry oni nit mesh This command is just shorthand for the calls pdeplot p e t pdeplot p e t zdata u If you want to have more control over your mesh plot usepdepl ot instead of pdemesh Plot the mesh for the geometry of the L shaped membrane p e t initmesh Ishapeg p e t refinemesh I shapeg p e t pdemesh p e t Now solve Poisson s equation Au 1 over the geometry defined by the L shaped membrane Use Dirichlet boundary conditions u 0 on Q and plot the result usassempde Ishapeb p e t 1 0 1 pdemesh p e t u pdeplot pdeco
73. 1 Q1 45Q2 5Q3 5Q4 C1 C2 C3 C4 1 29 1 Tutorial 1 30 Enter the set formula into the edit box at the top of the GUI Then enter the Boundary mode by pressing the dQ button or by selecting the Boundary Mode option from the Boundary menu The CSG model is now decomposed using the set formula and you get a rectangle with rounded corners as shown below xf PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help alol gt ac Pe A AN ALS Generic Soatar x 01 Y 06 Set formula ED 92 3 II cd 0 44 if 4 0 5L A IE a pi i i PE i pta i i i i i i i i i i 1 0 9 0 8 0 7 0 6 0 5 0 4 0 3 0 2 0 1 0 01 02 03 04 05 06 07 08 09 1 11 Info Push to define boundary cond s Exit Because of the intersection of the solid objects used in theinitial CSG model a number of subdomain borders remain They are drawn using gray lines If this is a model of e g a homogeneous plate you can remove them Select the Remove All Subdomain Borders option from the Boundary menu The subdomain borders are removed and the model of the plate is now complete Using the Graphical User Interface Suggested Modeling Method Although the PDE Toolbox offers you a great deal of flexibility in the ways that you can approach the problems and interact with the toolbox functions there is a Suggested method of choice
74. 139 Y 0 7782 f 1 1 f L L 1 L 1 L 1 14 12 1 0 8 06 04 0 2 0 0 2 0 4 06 0 8 1 12 14 Info Select a different plot or change mode to alter PDE mesh or boundary conditions Exit i The current density between two metallic conductors Heat Transfer The heat equation is a parabolic PDE T pC amp V KYT Q h T It describes the heat transfer process for plane and axi symmetric cases and uses the following parameters T ext e Density p e Heat capacity C e Coefficient of heat conduction k e Heat source Q e Convective heat transfer coefficient h e External temperature T ext 2 57 2 Examples The term h Teg T is a model of transversal heat transfer from the surroundings and it may be useful for modeling heat transfer in thin cooling plates etc For the steady state case the elliptic version of the heat equation V kVT O h T T is also available The boundary conditions can be of Dirichlet type where the temperature on the boundary is specified or of Neumann type where the heat flux n kV T is specified A generalized Neumann boundary condition can also be used The generalized Neumann boundary condition equation is n kV T qT g where q is the heat transfer coefficient Visualization of the temperature the temperature gradient and the heat flux kVT is available Plot options include isotherms and heat f
75. A and its residual is This residual has to be small in norm for 6 to be a good eigenvalue approximation The norm of the residual is rj hj 1 S the product of the last subdiagonal element of the Hessenberg matrix and the last element of its eigenvector It seldom happens that hj j gets particularly small but after sufficiently many steps j there are always some eigenvectors s with small last elements Note that the long vector Vj 41 is of unit norm It is not necessary to actually compute the eigenvector approximation y to get the norm of the residual we only need to examine the short vectors s and flag those with tiny last components as converged In a typical case n may be 2000 while j seldom exceeds 50 so all computations that involve only matrices and vectors of size j are much cheaper than those involving vectors of length n This eigenvalue computation and test for convergence is done every few steps j until all approximations to eigenvalues inside the interval Ib ub are flagged as converged When n is much larger than j this is done very often for smaller n more seldom When all eigenvalues inside the interval have converged or when j has reached a prescribed maximum the converged eigenvectors or more appropriately Schur vectors are computed and put in the front of the basis V After this the Arnoldi algorithm is restarted with a random vector if all approximations inside the interval are flagged as conv
76. ANd p e t poi mesh g USesnx ny 1 The triangular mesh is described by the mesh data p e andt Details on the mesh data representation can be found in the entry oni nit mesh For best performance with poi sol v the larger of nx and ny should be a power of 2 If y does not seem to describe a rectangle p is zero on return Try the demo command pdedemo8 The solution of Poisson s equation over a rectangular grid with boundary condition given by thefilesquareb4 is returned The solution time is compared to the usual FEM approach initmesh poisolv poisolv Purpose Synopsis Description See Also Reference Fast solution of Poisson s equation on a rectangular grid u poisolv b p e t f u poisolv b p e t f solves Poisson s equation with Dirichlet boundary conditions on a regular rectangular grid A combination of sine transforms and tridiagonal solutions is used for increased performance The boundary conditions b must specify Dirichlet conditions for all boundary points The mesh p e andt must be a regular rectangular grid Details on the mesh data representation can be found in the entry oni nit mesh f gives the right hand side of Poisson s equation Apart from roundoff errors the result should be the same as u assempde b p e t 1 0 f poicalc poi mesh Strang Gilbert Introduction to Applied Mathematics Wellesley Cambridge Press Cambridge MA 1986 pp 453 458 5 89 refinemesh Purpose Syn
77. DE itself through a dialog box that is accessed by pressing the button with the PDE icon or by selecting PDE Specification from the PDE pull down menu In the PDE mode you can also access the PDE Specification dialog box by double clicking on a subdomain That way different subdomains can have different PDE coefficient values This problem however consists of only one subdomain In the dialog box you can select thetype of PDE elliptic parabolic hyperbolic or eigenmodes and define the applicable coefficients depending on the PDE type This problem consists of an elliptic PDE defined by the equation V cVu au f with c 1 0 a 0 0 and f 10 0 G etting Started Finally create the triangular mesh that the PDE Toolbox uses in the Finite Element Method FEM tosolve the PDE Thetriangular mesh is created and displayed when pressing the button with the 44 icon or by selecting the Mesh menu option Initialize Mesh f you want a more accurate solution the mesh can be successively refined by pressing the button with the four triangle icon the Refine button or by selecting the Refine Mesh option from the Mesh menu Using theJ iggle Mesh option the mesh can be jiggled to improve the triangle quality Parameters for controlling the jiggling of the mesh the refinement method and other mesh generation parameters can be found in a dialog box that is opened by selecting Parameters from the Mesh menu You can undo any chan
78. E and applicable coefficients Exit Visualization of the von Mises effective stress and the displacements using deformed mesh Structural Mechanics Plane Strain A deformation state where there are no displacements in the z direction and the displacements in the x and y directions are functions of x and y but not z is called planestrain You can solve plane strain problems with the PDE Toolbox by selecting the application mode Structural Mechanics Plane 2 41 2 Examples Strain The stress strain relation is only slightly different from the plane stress case and thesame set of material parameters is used The application interfaces are identical for the two structural mechanics modes The places where the plane strain equations differ from the plane stress equations are e Theu parameter in the c tensor is defined as v 1 2v e The von Mises effective stress is computed as u 2G NV v 1 GIG N 1 Plane strain problems are less common than plane stress problems An example is a slice of an underground tunnel that lies along the z axis It deforms in essentially plane strain conditions 2 42 Application Modes Electrostatics Applications involving electrostatics include high voltage apparatus electronic devices and capacitors The statics implies that the time rate of change is slow and that wavelengths are very large compared to the size of the domain of interest In elec
79. E solution vector u and if applicable the computed eigenvalues to the main workspace Parameters Solve Parameters dialog box for elliptic PDEs 3 25 3 The Graphical User Interface 3 26 Parameters opens a dialog box where you can enter the solve parameters The set of solve parameters differs depending on the type of PDE e Elliptic PDEs By default no specific solve parameters are used and the elliptic PDEs are solved using the basic elliptic solver ass e mp de Optionally the adaptive mesh generator and solver adapt mesh can be used For the adaptive mode the following parameters are available Adaptive mode Toggle the adaptive mode on off Maximum number of triangles The maximum number of new triangles allowed can beset to nf A default value is calculated based on the current mesh Maximum number of refinements The maximum number of successive refinements attempted Triangle selection method Two triangle selection methods are included in the PDE Toolbox You can also supply your own function Worst triangles This method picks all triangles that are worsethan a fraction of the value of the worst triangle default 0 5 SeetheReference entry pdeadworst for more details Relative tolerance This method picks triangles using a relative toler ance criterion default 1E 3 SeetheReferenceentry pdeadgsc for more details User defined function Enter the name of a user defined triangle sel
80. M fileiscirclef that produces a point load Y ou can find thefull exampleinpdedemo7 andintheReferenceentry onassempde The routines for adaptive mesh generation and solution are powerful but can lead to dense meshes and thus long computation times Setting the Nge n parameter to one limits you to a single refinement step This step can then be repeated to show the progress of the refinement TheMaxt parameter helps you stop before the adaptive solver generates too many triangles An example of a hand written triangle selection function iscirclepick usedinpdedemo7 Remember that you always need a decomposed geometry with adapt mesh Deformed meshes are easily plotted by adding offsets to the Point matrix p Assuming two variables stored in the solution vector u np size p 2 pdemesh ptscale u l np u nptl nptnp e t The time evolution of eigenmodes is obtained by e g ul u mode cos sqrt I mode tlist hyperbolic for positive eigenvalues in hyperbolic problems or ul u mode exp I mode tlist parabolic in parabolic problems This makes nice animations perhaps together with deformed mesh plots 1 41 1 Tutorial 1 42 Examples Examples of Elliptic Problems 2 2 Poisson s Equation on Unit DiSk o o oooooooomoomoo 2 2 A Scattering Problem o ooocoocorcoo 2 6 A Minimal Surface Problem 00 02 eee eee eee 2 10 Domain Decomposition oooooocoo
81. N 2N 1 2 Yes Gikl 2j j 1 k 4 Gjki 2j 3j 4i 21 k 5 4N No Gjki AN j 1 4i 21 k 6 An example of how c is stored inc can be found below assempde Examples Example 1 Solve the equation Au 1 on the geometry defined by the L shaped membrane Use Dirichlet boundary conditions u 0 on 0Q Finally plot the solution p e t initmesh Ishapeg Hmax 0 2 p e t refinemesh Ishapeg p e t usassempde Ishapeb p e t 1 0 1 pdesurf p t u Example 2 Consider Poisson s equation on the unit circle with unit point source at the origin The exact solution u L log r is known for this problem We define the function f circlef p t u time for computing the right hand side circlef returns zero for all triangles except for the one located at the origin for that triangle it returns 1 a wherea isthetrianglearea pdedemo7 performs a full demonstration of the problem with adaptive solution 5 23 assempde 5 24 Example 3 We study how the matrices a and also d are stored in the MATLAB matrix a for system case N 3 Na 1 The bullet symbol below means that the matrix is symmetric na 6 ne 1 N 3 Ne 2 a 1 Na 3 a 1 Na 9 a 1 a 2 a 4 a 3 a 5 e e a 6 a 1 0 a 2 0 a 3 a 1 a 4 a 7 a 2 a 5 a 8 a 3 a 6 a 9
82. Name Space Matrix The Name Space matrixns relates the columns in gd to variable names in sf Each column in ns contains a sequence of characters padded with spaces Each such character column assigns a name tothe corresponding geometric object in gd This way we can refer toa specific object in gd in the set formula sf Decomposed Geometry Matrix The Decomposed Geometry matrix dI contains a representation of the decomposed geometry in terms of disjointed minimal regions that have been constructed by thedecsg algorithm Each edge segment of the minimal regions corresponds to a column in dl We refer to edge segments between minimal regions as border segments and outer boundaries as boundary segments In each such column rows two and three contain the starting and ending x coordinate and rows four and five the corresponding y coordinate Rows six and seven contain left and right minimal region labels with respect to the direction induced by the start and end points counter clockwise direction on decsg circle and ellipse segments There are three types of possible edge segments in a minimal region e For circle edge segments row one is 1 Rows eight and nine contain the coordinates of the center of the circle Row 10 contains the radius For line edge segments row one is 2 For ellipse edge segments row one is 4 Rows eight and nine contain the coordinates of the center of the ellipse Rows 10 and 11 contain the semiaxes of the ellips
83. Neumann boundary condition requires the value of the normal component of the current density n o V V to be known It is also possible to specify a generalized Neumann condition defined by n o V V qV g where q can beinterpreted as a film conductance for thin plates The electric potential V the electric field E and the current density J areall available for plotting Interesting quantities to visualize are the current lines the vector field of J and the equipotential lines of V The equipotential lines are orthogonal to the current lines when o is isotropic Example Two circular metallic conductors are placed on a plane thin conductor likea blotting paper wetted by brine The equipotentials can be traced by a voltmeter with a simple probe and the current lines can be traced by strongly colored ions such as permanganate ions The physical model for this problem consists of the Laplace equation V oVV 0 for the electric potential V and the boundary conditions e V 1 on the left circular conductor e Y 1ontheright circular conductor e Thenatural Neumann boundary condition Y 0 on the outer boundaries n The conductivity o 1 constant 2 55 2 Examples Using the Graphical User Interface Inthepdetool GUI select the application mode Conductive Media DC and draw the blotting paper as a rectangle with corners in 1 2 0 6 1 2 0 6 1 2 0 6 and 1 2 0 6 Add two circles with radius 0 3 represe
84. Partial Differential Equation Toolbox Fee perp AE For Use with MATLAB Computer Solutions Europe AB Computation Visualization i Programming User s Guide X OD o E How to Contact The MathWorks 508 647 7000 Phone 508 647 7001 Fax The MathWorks Inc Mail 24 Prime Park Way Natick MA 01760 1500 http www mathworks com Web ftp mathworks com Anonymous FTP server comp soft sys matlab Newsgroup support mathworks com Technical support suggest mathworks com Product enhancement suggestions bugs mathworks com Bug reports doc mathworks com Documentation error reports subscribe mathworks com Subscribing user registration service mathworks com Order status license renewals passcodes info mathworks com Sales pricing and general information Partial Differential Equation Toolbox User s Guide COPYRIGHT 1984 1997 by The MathWorks Inc All Rights Reserved The software described in this document is furnished under a license agreement The software may be used or copied only under the terms of the license agreement No part of this manual may be photocopied or repro duced in any form without prior written consent from The MathWorks Inc U S GOVERNMENT If Licensee is acquiring the Programs on behalf of any unit or agency of the U S Government the following shall apply a For units of the Department of Defense the Government shall have only the rights specified in the license under which the comm
85. Synopsis Description See Also Indices of points in canonical ordering for rectangular grid nl n2 hl h2 i1 c ii cc poiindex p e t sd nl n2 h1 h2 i c ii cc poiindex p e t sd identifies a given gridp e t in the subdomain sd as an evenly spaced rectangular grid If the grid is not rectangular n1 isO on return Otherwisen1 andn2 arethe number of points in the first and second directions h1 andh2 arethespacings i andi are of length n1 2 n2 2 and contain indices of interior points contains indices of the original mesh whereas i i contains indices of the canonical ordering c and cc are of length n1 n2 n1 2 n2 2 and contain indices of border points i i andcc areincreasing In the canonical ordering points are numbered from left to right and then from bottom to top Thus if n1 3 andn2 5 thenii 5 8 11 and cc 1 23 467 9 10 12 13 14 15 poisolv poiasma 5 87 poimesh Purpose Synopsis Description Examples See Also 5 88 Make regular mesh on a rectangular geometry p e t poi mesh g nx ny p e t poi mesh g n p e t poi mesh g p e t poi mesh g nx ny Constructs a regular mesh on the rectangular geometry specified by g by dividing the x edge intonx pieces and the y edge intony pieces and placing nx 1 ny 1 points at the intersections The x edge is the one that makes the smallest angle with the x axis p e t poi mesh g n Usesnx ny n
86. The triangles form a mesh and each vertex is called a node Y ou are in the situation of an architect designing a dome He has to strike a balance between the ideal rounded forms of the original sketch and the limitations of his simple building blocks triangles or quadrilaterals If the result does not look close enough to a perfect dome the architect can always improve his work using smaller blocks Next you say that your solution should besimple on each triangle Polynomials are a good choice they are easy to evaluate and have good approximation properties on small domains You can ask that the solutions in neighboring triangles connect to each other continuously across the edges Y ou can still decide how complicated the polynomials can be J ust like an architect you want them as simple as possible Constants are the simplest choice but you cannot match values on neighboring triangles Linear functions come next Basics of The Finite Element M ethod This is like using flat tiles to build a waterproof dome which is perfectly possible DAD i Heh ABS TAN YN ANY Ns AN sx A triangular mesh left and a continuous piecewise linear function on that mesh Now you use the basic elliptic equation V cVu au f ing If u is the piecewise linear approximation to u it is not clear what the second derivative term means Inside each triangle Vu is a constant because u is flat and thus the secon
87. Then define the boundary conditions by pressing the dQ button and then double click on the boundaries to define the boundary conditions On the right hand side boundary you have the generalized Neumann conditions and you enter them as constants g 0 and q 3 4 Initialize the mesh and refine it once by pressing the A and refine buttons or by selecting the corresponding options from the Mesh menu Also define the eigenvalue PDE problem by opening the PDE Specification dialog box and selecting the Eigenmodes option The general eigenvalue PDE is described by V cVu au Adu so for this problem you use the default values c 1 a 0 and d 1 Also in the Solve Parameters dialog box enter the eigenvalue range as the MATLAB vector I nf 10 Finally press the button to compute the solution By default the first eigenfunction is plotted Y ou can plot the other eigenfunctions by selecting the corresponding eigenvalue from a pop up menu in the Plot Selection dialog box The pop up menu contains all the eigenvalues found in the specified range You can also export the eigenfunctions and eigenvalues to the MATLAB main workspace by using the Export Solution option from the Solve menu Using Command Line Functions The geometry description file and boundary condition file for this problem are calledsquareg mandsquareb2 m respectively Use the following sequence of commands to find the eigenvalues in the specified range and the corr
88. There Any Applications Already Implemented 1 7 Can Extend the Functionality of the Toolbox 1 7 How Can Solve 3 D Problems by 2 D Models 1 8 Getting Started 0 0 eee 1 9 Basics of The Finite Element Method 1 18 Using the Graphical User Interface 1 23 The PDE Toolbox Graphical User Interface 1 23 The M NUS session rr a ow a a ees eka 1 24 THO TOO DAR Sent feck Bw cows Be aa ah OE we ae Hae ecw OS 1 25 The GULModes vivian tea tie Sawa WaT Wace wR no aa 1 26 The CSG Model and the Set Formula o 1 27 Creating Rounded Corners oooccccccc nee 1 28 Suggested Modeling Method 0 cece eee ee eens 1 31 Object Selection Methods 0 cece ee eee 1 35 Display Additional Information o oocooooooooo 1 35 Entering Parameter Values as MATLAB Expressions 1 36 Using PDE Toolbox version 1 0 Model M files 1 36 Using Command Line Functi0nS 1 37 Data Structures and Utility Functions 1 37 Constructive Solid Geometry Model 1 38 Decomposed Geometry o ooocccoccoc tee ee 1 39 Boundary Conditions 0 0c eee eee 1 39 Equation Coefficients 0 cece ees 1 39 MCSA score a hod Ae wh wm Kose wc awh apa Re ata T i e d a A i 1 39 SOLO 28 ade ied Pale a dala aes ae Dada awe eds 1 40 Post Processing and Pres
89. a 1 19 4 4 Grid Spacing command 3 10 H heat distribution in radioactive rod 2 21 heat equation 2 16 2 57 4 13 heat transfer 2 57 Helmholtz s equation 2 6 2 52 Help menu 3 37 hyperbolic 5 36 hyperbolic equation 1 3 4 16 hyperbolic problems 2 23 l initmesh 5 38 jigglemesh 5 42 L Laplace equation 2 44 2 55 L shaped membrane 2 27 2 31 M magnetostatics 2 46 mass matrix 4 6 4 16 Maxwell s equations 2 43 2 46 2 51 menus 1 24 3 3 mesh data 1 39 5 38 Mesh menu 3 22 mesh mode 1 27 mesh parameters 3 23 method of lines 4 14 minimal region 1 39 5 32 minimal surface problem 2 10 mixed boundary condition 4 3 4 11 Model M file 1 36 5 64 5 81 N Name Space matrix 1 38 5 32 Neumann boundary condition 4 3 4 11 New command 3 3 node data 5 62 5 70 nonlinear equation 4 21 nonlinear problem 2 10 nonlinear solver 1 21 0 Open command 3 4 Options menu 3 9 Index P parabolic 5 43 parabolic equation 1 3 4 13 parabolic problems 2 16 Paste command 3 8 PDE coefficients 3 19 5 20 5 21 PDE menu 3 18 PDE mode 1 27 PDE Specification 3 19 pdeadgsc 5 45 pdeadworst 5 46 pdearcl 5 47 pdebound 5 48 pdecgrad 5 50 pdecirc 5 51 pdecont 5 52 pdeeig 5 53 pdeellip 5 55 pdegeom 5 57 pdegplot 5 60 pdegrad 5 61 pdeintrp 5 62 pdejmps 5 63 pdemdl cv 5 64 pdemesh 5 65 pdenonlin 5 66 pdeplot 5 69 pdepol y 5 71 pdeprtni 5 72 pderect 5 73 pdesde 5 74 pdesdp 5 74 pdesdt 5 74
90. ains Calculate structural mechanics tensor functions Triangle geometry data 5 5 5 Reference Function Purpose pdetriq Triangle quality measure poi asma Boundary point matrix contributions for fast solvers of Poisson s equation poicalc Fast solver for Poisson s equation on a rectangular grid poiindex Indices of points in canonical ordering for rectangular grid sptarn Solve generalized sparse eigenvalue problem tri2grid Interpolate from PDE triangular mesh to rectangular grid Commands G rouped by Function User Defined Algorithms Function Purpose pdebound Boundary M file pdegeom Geometry M file Demonstration Programs Function Purpose pdedemol Exact solution of Poisson s equation on unit disk pdedemo2 Solve Helmholtz s equation and study the reflected waves pdedemo3 Solve a minimal surface problem pdedemo4 Solve PDE problem using subdomain decomposition pdedemo5 Solve a parabolic PDE the heat equation pdedemo6 Solve a hyperbolic PDE the wave equation pdedemo7 Adaptive solution with point source pdedemo8 Solve Poisson s equation on rectangular grid 5 7 adaptmesh Purpose Synopsis Description 5 8 Adaptive mesh generation and PDE solution u p e t adaptmesh g b c a f u p t adaptmesh g b c a f PropertyName PropertyValue u p e t adaptmesh g b c a f PropertyName PropertyValue performs adaptive mesh gene
91. alyzes the CSG model constructs a set of disjoint minimal regions bounded by boundary segments and border segments and optionally evaluates a set formula in terms of the objects in the CSG model We often refer to the set of minimal regions as the decomposed geometry The decomposed geometry makes it possible for other toolbox functions to understand the geometry you specify For plotting purposes a second set of minimal regions with a connected boundary is constructed The graphical user interfacepdetool usesdecsg for many purposes Each time a new solid object is drawn or changed pdetool callsdecsg in order to beable to draw the solid objects and minimal regions correctly The Delaunay triangulation algorithm i ni t mesh also uses the output of de c s y to generate an initial mesh dl decsg gd sf ns decomposes the CSG model gd into the decomposed geometry d The CSG model is represented by the Geometry Description matrix and the decomposed geometry is represented by the Decomposed Geometry matrix decs g returns the minimal regions that evaluate to true for the set formulas f The Name Space matrix ns is a text matrix that relates the columns in gd to variable names in sf dl decsg gd returnsall minimal regions The same as letting sf correspond to the union of all objects in gd dl bt decsg gd and dl bt decsg gd sf ns additionally return a Boolean table that relates the original solid objects to the minimal region
92. an use any unique name as long as it contains no blanks In Draw mode you can alter the names and the geometries of the objects by double clicking on them This opens a dialog box where you can 1 27 1 Tutorial 1 28 edit the name and the geometry The following figure shows an object dialog box for a circle ir Object Dialog Object type Circle X center 0 5 Y center 0 3 Radius 15 Name C2 OK Cancel ol You can use the name of the object to refer to the corresponding set of points in a set formula The operators and are used to form the set of points Q in the plane over which the differential equation is solved The operators the set union operator and the set intersection operator have the same precedence The operator the set difference operator has higher precedence The precedence can be controlled by using parentheses The resulting geometrical model Q is the set of points for which the set formula evaluates to true By default it is the union of all solid objects We often refer to the area Q as the decomposed geometry Creating Rounded Corners As an example of how to use the set formula let s model a plate with rounded corners fillets Start the GUI and turn on the grid and the snap to grid feature using the Options menu Also change the grid spacing to 1 5 0 1 1 5 for the x axis and 1 0 1 1 for the y axis Select Rectangle square from the Dr
93. arly visible as it is drawn using a darker shade of gray The object that you just drew the circle has a black border indicating that it is selected A selected object can be moved resized copied and deleted Y ou can select more than one object by Shift clicking on the objects that you want to select Also a Select All option is available from the E dit menu x PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help 1 F O eS gt aa PDE Z 4A A X Genetic Soatar al X 15 Y 0 8 Set formula R1 C1 R2 C2 cl R1 1 i i i i i 15 1 0 5 a 0 5 1 15 Info Draw and edit 2 D geometry by using the Draw and Edit menu options Exit g y by g pl Finally add two more objects a rectangle R2 and a circle C2 The desired CSG model is formed by subtracting the circle C2 from the union of the other three objects You dothis by editing the set formula that by default is the union of all objects C1H14R24C2 You can type any other valid set formula into Set formula edit field Click in the edit field and use the keyboard to change the set formula to R1 C1 R2 C2 1 11 1 Tutorial If you want you can save this CSG model as an M file Use the Save As option from the File menu and enter a filename of your choice It s good practice to continue to save your model at regular intervals using Save All the addi
94. ate The rate at which the mesh size increases away from small parts of the geometry The value must be between 1 and 2 The default value is 1 3 i e the mesh size increases by 30 e J iggle mesh Toggles automatic jiggling of the initial mesh on off 3 23 3 The Graphical User Interface 3 24 The parameters used by the mesh jiggling algorithm i gg emesh are e J iggle mode Select a jiggle mode from a pop up menu Available modes are on optimize minimum andoptimize mean on jiggles the mesh once Using the jiggle modeoptimize mi ni mum thejiggling process is repeated until the minimum triangle quality stops increasing or until the iteration limit is reached The same applies for theopti mize mean option but it tries to increase the mean triangle quality Number of jiggle iterations Iteration limit for theopti mize mini mum and optimize mean modes Default 20 Finally for the mesh refinement algorithm efi nemesh the Refinement method can beregu ar or ongest The default refinement method is regul ar which results in a uniform mesh The refinement method ongest always refines the longest edge on each triangle PDE Toolbox M enus Solve Menu Solve PDE Solve the partial differential equation for the current CSG model and triangular mesh and plot the solution the automatic solution plot can be disabled Parameters Open dialog box for entry of PDE solve parameters Export Solution Export the PD
95. aw menu or press the button with the rectangle icon Then draw a rectangle with a width of 2 and a height of 1 using the mouse starting at 1 0 5 To get the round corners add circles onein each corner The circles should have a radius of 0 2 and centers at a distance that is 0 2 units from the left right and lower upper rectangle boundaries 0 8 0 3 0 8 0 3 0 8 0 3 and 0 8 0 3 To draw several circles double click on the button for drawing ellipses circles centered Then draw the circles using the right mouse button or Control click starting at the circle centers Finally at each of the rectangle corners draw four small squares with a side of 0 1 Using the Graphical User Interface The figure below shows the complete drawing PDE Toolbox Untitled ile Edit Options Draw Boundary PDE Mesh Solve Plot Window Help Set formula R14014024C3 C44901 490248034904 0 5 sez cl og 0 1 R1 4 SQ3 c3 cd sq 4 0 54 4 4 0 7 f E A i i i i i i i i i i i i i i i i i i i i i 1 0 9 0 8 0 7 0 6 0 5 0 4 0 3 0 2 0 1 0 01 02 03 04 05 06 07 08 09 1 11 Info Click at corner and drag to create rectangle square Exit Now you have to edit the set formula To get the rounded corners subtract the small squares from the rectangle and then add the circles As a set formula this is expressed as R
96. ay name clashes between variables used in the function and in the main workspace are avoided The name of the file must coincide with the model name The beginning of the file always looks similar to the code fragment below function pdemode pdeinit pde_fig gcf ax gca pdetool appl_cb 1 setuprop pde_fig currparam str2mat 1 0 0 0 10 0 1 0 pdetool snapon set ax AspectRatio 1 5 1 set ax XLim 1 5 1 5 set ax YLim 1 1 set ax XTickMode auto set ax YTickMode auto grid on Thepdei nit command starts uppdetool Ifpdetoo has already been started the current model is cleared The following commands set up the scaling and tick marks of the axis of pdetoo l and other user parameters Then a sequence of drawing commands is issued The commands that can be used arenamed pdecirc pdeellip pdepol y andpder ect The command sequence below creates the L shaped membrane as the union of three squares The solid objects are given names SQ1 SQ2 SQ3 etc Geometry description pderect 1 0 0 1 SQ1 pderect 0 1 0 1 SQ2 pderect 0 1 1 0 SQ3 5 81 pdetool We do not intend to fully document the format of the M odel M file It can be used to change the geometry of the drawn objects sincethepdecirc pdeellip pdepoly andpderect commands are documented See Also initmesh assempde parabolic wave pdeeig pdesurf pdecont 5 82 pdetrg Purpose Synopsis Descr
97. b p e t c a f time sdl pay E a f u0 time e t c B ud assempde b p e t c a f u0 e t c ud assempde b p e t c a f ti me Q G H R assempde b p e t c a f Q G H R assempde b p e t c a f Uu0 Q G H R assempde b p e t c a f u0 ti me Q G H b p e t Q G b p e t b p e t R assempde c a f u0 time sdl H R assempde C a f time K M F Q G H R assempde c a f time sdi u assempde K M F Q G H R K1 Fl assempde K M F Q G H R K1 F1 B ud assempde K 1 1 1 assempde isthebasicfunction of the PDE Toolbox It assembles a PDE problem by using the FEM formulation described in Chapter 4 The Finite Element Method The command assempde assembles the scalar PDE problem V cVu au f on Q or the system PDE problem V c 8 Vu au f on Q The command can optionally produce a solution tothe PDE problem assempde For the scalar case the solution vector u is represented as a column vector of solution values at the corresponding node points from p For a system of dimension N with np node points the first np values of u describe the first component of u the following np values of u describe the second component of u and soon Thus the components of u are placed in the vector u as N blocks of node point values u assempde b p e t c a f assembles and solves the PDE problem by eliminating the Dirichlet boundary conditions from the system of linear equations
98. box where you can define the PDE coefficientsc a andf Inthis simple case they are all constants c 1 f 1 anda 0 Press the 4 button or select Initialize Mesh from the Mesh menu This initializes and displays a triangular mesh Press the Refine button or select Refine Mesh from the Mesh menu This causes a refinement of the initial mesh and the new mesh is displayed Tosolvethe system just press the button The toolbox assembles the PDE problem and solves the linear system It also provides a plot of the solution Using the Plot Selection dialog box you can select different types of solution plots xj PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help S gt a P Al A AJA Generic Scalar AEREE Y 06 Set formula k Color u Into Select a different plot or change mode to alter PDE mesh or boundary conditions To compare the numerical solution to the exact solution select the user entry in the Property pop up menu for Color in the Plot Selection dialog box Then input the MATLAB expression u 1 x 2 y 2 4 in the user entry edit field You obtain a plot of the absolute error in the solution 2 Examples 2 4 You can also compare the numerical solution to the exact solution by entering some simple command line oriented commands see the next section It is easy to export the mesh data an
99. c length values p is a monotone row vector of parameter values and x y is a matrix with two rows giving the corresponding points on the curve The first point of the curve is given the arc length values 0 and the last point the values 1 On return pp contains parameter values corresponding to the arc length values specified ins The arc length valuess s0 ands1 can be an affine transformation of the arc length See the example car dg in the entry pdegeom pdegeom pdebound Purpose Synopsis Description Boundary M file q 9 h r pdebound p e u ti me The Boundary M file specifies the boundary conditions of a PDE problem The most general form of boundary conditions that we can handleis hu r CQ Vu qu g p By the notation 7 c Vu we mean the N by 1 matrix with i 1 component N lo lo lo lo y cos amp c 113 cos 0 c 1255 sin amp c a sin O y Ui Jal where the outward normal vector of the boundary cos a sin a There are M Dirichlet conditions and the h matrix is M by N M gt 0 The generalized Neumann condition contains a source k u where the Lagrange multipliers u is computed such that the Dirichlet conditions become satisfied The data that you specify is q g h andr For M 0 we say that we have a generalized Neumann boundary condition for M N a Dirichlet boundary condition and for 0 lt M lt N a mixed boundary condition The Boundary M
100. case you must use your knowledge of the original physical problem to interpret the results from the computation sptarn pdeellip Purpose Synopsis Description Examples See Also Draw ellipse pdeellip xc yc a b phi pdeellip xc yc a b phi label pdeellip xc yc a b phi draws an ellipse with center in xc yc and semiaxes a andb Therotation of the ellipse in radians is given by phi If the pdetool GUI is not active it is automatically started and the ellipse is drawn in an empty geometry model The optional argument abel assigns a nametotheellipse otherwisea default name is chosen The state of the Geometry Description matrix inside pdetool is updated to include the ellipse You can export the Geometry Description matrix from pdetool by selecting the Export Geometry Description option from the Draw menu Theformat of the Geometry Description matrix is described in the entry on decsg The command below starts pdetoo and draws an ellipse pdeellip 0 0 1 0 3 pi 4 pdecirc pdepoly pderect pdetool 5 55 pdeent Purpose Indices of triangles neighboring a given set of triangles Synopsis ntl pdeent t tl Description Given triangle datat anda list of triangleindicest nt contains indices of the triangles int and their immediate neighbors i e those whose intersection witht is nonempty See Also refinemesh 5 56 pdegeom Purpose Synopsis Description Geometry
101. center Then put the cursor at the desired corner and click and drag using the left mouse button to create a rectangle with the desired side lengths Use the right mouse button to create a square Notice how the snap to grid feature forces the rectangle to line up with the grid When you release the mouse the CSG model is updated and redrawn At this stage all you have is a rectangle It is assigned the label R1 If you want to move or resize the rectangle you can easily do so Click and drag an object to move it and double click on an object to open a dialog box where you can enter exact location coordinates From the dialog box you can also alter the label If you are not satisfied and want to restart you can delete the rectangle by pressing the Delete key or by selecting Clear from the Edit menu Next draw a circle by clicking on the button with the ellipse icon with the sign and then click and drag in a similar way using the right mouse button starting at the circle center xf PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help 0 2b Ri 7 1 i i i i i 15 1 0 5 0 Os 1 15 Info Draw and edit 2 D geometry by using the Draw and Edit menu options Exit i 1 10 G etting Started The resulting CSG model is the union of the rectangle R1 and the circle C1 described by set algebra as R1 C1 The area where the two objects overlap is cle
102. ch covers all aspects of the PDE solution process You start it by typing pdetool at the MATLAB command line It may take a while the first time you launch pdetool during a MATLAB session The figure below shows thepdetoo GUI as it looks when you have started it Y PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help a gt gt 30 poe A Ax il A Generic Scalar Xx 00 Y 0 0 Set formula 1 T T T ash 4 0 6 H 4 0 41 4 02 4 ob H 0 24 4 0 44 4 0 6L El 0 8L el 1 1 1 1 1 1 15 1 0 5 0 Os 1 15 f Info Click and drag at perimeter to create ellipse Exit A A At the top the GUI has a pull down menu bar that you use to control the modeling It conforms to common pull down menu standards Menu items followed by a right arrow lead toa submenu Menu items followed by an ellipsis lead to a dialog box Stand alone menu items lead to direct action Below the menu bar a toolbar with icon buttons provide quick and easy access to some of the most important functions 1 23 1 Tutorial Tothe right of the toolbar is a pop up menu that indicates the current application mode You can also use it to change the application mode The upper right part of the GUI also provides the x and y coordinates of the current cursor position It is updated when you move the cursor inside the main axes area in the middle of
103. ch triangle into four similar triangles by creating new corners at the midsides adjusting for curved boundaries You can assess the accuracy of the numerical solution by comparing results from a sequence of successively refined meshes If the solution is smooth enough more accurate results may be obtained by extrapolation The solutions of the toolbox equation often have geometric features like localized strong gradients An example of engineering importance in elasticity is thestress concentration occurring at re entrant corners such as the MATLAB favorite the L shaped membrane Then it is more economical to refine the mesh selectively i e only where it is needed When the selection is based on estimates of errors in the computed solutions a posteriori estimates we speak of adaptive mesh refinement See the Reference entry on adapt mesh for an example of the computational savings where global refinement needs more than 6000 elements to compete with an adaptively refined mesh of 500 elements The adaptive refinement generates a sequence of solutions on successively finer meshes at each stage selecting and refining those elements that are judged to contribute most to the error The process is terminated when the maximum number of elements is exceeded or when each triangle contributes less than a preset tolerance Y ou need to provide an initial mesh and choose selection and termination criteria parameters The initial mesh can be produced by
104. click click and drag to create a square Rectangle square centered Draw a rectangle square starting at the center Using the left mouse button click and drag to create a rectangle Using the right mouse button or Control click click and drag to create a square Ellipse circle Draw an ellipse circle starting at the perimeter Using the left mouse button click and drag to create an ellipse Using the right mouse button or Control click dick and drag to create a circle 3 13 3 The Graphical User Interface Ellipse circle centered Polygon Rotate Export Geometry Description Set Formula Labels Draw an ellipse circle starting at the center Using the left mouse button click and drag to create an ellipse Using the right mouse button or Control click click and drag to create a circle Draw a polygon You can close the polygon by pressing the right mouse button Clicking at the starting vertex also closes the polygon Rotate selected objects Export the Geometry Description matrix gd the set formula string sf and the Name Space matrixns labels to the main workspace Rotate Rotation degrees E Use center of mass Cancel all Rotate opens a dialog box where you can enter the angle of rotation in degrees The selected objects are then rotated by the number of degrees that you specify The rotation is done counter clockwise for positive rotatio
105. command 3 11 B Boolean table 5 30 border segment 5 32 5 57 boundary condition 3 16 4 3 4 11 5 14 5 80 Boundary Condition matrix 1 39 5 14 Boundary menu 3 15 Boundary M file 1 39 5 48 boundary mode 1 26 boundary segment 5 32 5 57 C circle solid 5 31 Coefficient matrix 1 39 5 21 Coefficient M file 1 39 5 21 command line functions using 1 37 2 4 2 11 2 19 2 25 2 28 2 33 conductive media DC 2 55 Constructive Solid Geometry model 1 27 1 38 5 30 5 78 CSG model 1 27 1 38 5 30 5 78 csgchk 5 28 csgdel 5 29 cylindrical problem 2 21 D decomposed geometry 1 28 1 39 5 30 Decomposed Geometry matrix 1 39 5 32 decsg 5 30 Dirichlet boundary condition 4 3 4 11 discrete sine transform 4 29 5 35 domain decomposition 2 12 Draw menu 3 13 draw mode 1 26 dst 5 35 E E dge matrix 1 39 5 38 Edit menu 3 7 eigenmodes 2 32 eigenvalue equation 4 17 eigenvalue problems 1 21 2 27 5 53 eigenvector matrix 5 53 electrostatics 2 43 ellipse solid 5 32 elliptic equation 1 3 1 19 4 3 elliptic problems 2 2 elliptic system 4 10 energy norm 4 6 Index 1 2 F fast solver 4 29 5 89 FEM 1 18 4 2 File menu 3 3 Finite Element Method 1 18 4 1 G Gauss N ewton method 4 22 Geometry Description matrix 1 38 5 31 Geometry M file 1 39 5 57 graphical user interface using 1 9 1 23 2 2 2 8 2 10 2 17 2 22 2 23 2 27 2 33 2 39 2 44 2 48 2 53 2 56 2 59 Green s formul
106. ct The condition at the outer computational boundary should be chosen to allow waves to pass without reflection Such conditions are usually called nonreflecting and weusetheclassical Sommerfeld radiation condition As approaches infinity r approximately satisfies the one way wave equation or Vr 0 which allows waves moving in the positive direction only 8 is the radial distance from the object With the time harmonic solution this turns into the generalized Neumann boundary condition E E Vr ikr For simplicity reasons let s make the outward normal of the computational domain approximate the outward direction 2 Examples Using the Graphical User Interface You can now usepdetoo tosolvethis scattering problem Using the Generic Scalar mode start by drawing the 2 D geometry of the problem Let the illuminated object be a square SQ1 with a side of 0 1 units and center in 0 8 0 5 and rotated 45 degrees and let the computational domain be a circle C1 with a radius of 0 45 units and the same center location The Constructive Solid Geometry CSG model is then given by C1 5Q1 For the outer boundary the circle perimeter the boundary condition is a generalized Neumann condition with q ik The wave number k 60 which corresponds to a wavelength of about 0 1 units so enter 601 as a constant q and0 asa constant g For the square object s boundary you have a Dirichlet boundary condition ik
107. cted plot options If you haven t defined your PDE pdet ool solves the default PDE which is Poisson s equation A u 10 This corresponds to the generic elliptic PDE with c 1 a 0 and f 10 For the different application modes different default PDE settings apply 1 34 Using the Graphical User Interface Object Selection Methods Throughout the GUI similar principles apply for selecting objects such as solid objects subdomains and boundaries e Toselect a single object click on it using the left mouse button e To select several objects and to deselect objects Shift click or click using the middle mouse button on the desired objects e Clicking in the intersection of several objects selects all the intersecting objects e To open an associated dialog box double click on an object If the object is not selected it is selected before opening the dialog box e n Draw mode and PDE mode clicking outside of objects deselects all objects e To select all objects use the Select All option from the Edit menu e When defining boundary conditions and the PDE via the menu items from the Boundary and PDE menus and no boundaries or subdomains are selected the entered values applies to all boundaries and subdomains by default Display Additional Information In Mesh mode you can use the mouse to display the node number and the triangle number at the position where you click Press the left mouse button to di
108. d order term vanishes At the edges of the triangles cVu isin general discontinuous and a further derivative makes no sense What you are looking for is the best approximation of u in the class of continuous piecewise polynomials Therefore you test the equation for u against all possible functions v of that class Testing means formally to multiply the residual against any function and then integrate i e determine u such that fa EY cVu au f vdx 0 for all possible v The functions v are usually called test functions Partial integration Green s formula yields that u should satisfy fa cVu v au vdx oon cVu vds faa Vv 1 19 1 Tutorial where dQ is the boundary of Q and a is the outward pointing normal on dQ Note that the integrals of this formulation are well defined even if u and v are piecewise linear functions Boundary conditions are included in the following way If u is known at some boundary points Dirichlet boundary conditions we restrict the test functions tov O at those points and require u to attain the desired value at that point At all the other points we ask for Neumann boundary conditions i e cVu qu g The FEM formulation reads Find u such that fa eva Vv au Vdx lee qu vds faa ha 2vds Yv where 0Q is the part of the boundary with Neumann conditions The test functions v must be zero on 0Q 0Q Any continuous piecewise linear u is represen
109. d the solution to MATLAB s main workspace by using the Export options from the Mesh and Solve menus To refine the mesh and solve the PDE successively simply press the refine and push buttons until the desired accuracy is achieved Another possibility is tousethe adaptive solver Using Command Line Functions First you must createa MATLAB function that parameterizes the 2 D geometry inthis casea unit circle The M fileci rcl eg m returns the coordinates of points on the unit circle s boundary The file conforms to the file format described in the pdegeom entry in the Reference chapter You can display the M file by typingtype circleg Also you need a function that describes the boundary condition This is a Dirichlet boundary condition where u 0 on the boundary The M file circleb1 m provides the boundary condition The file conforms to the file format described in thepdebound entry in Chapter 5 Reference You can display the M file by typingt ype circlebl Now you can start working from MATLAB s command line p e t initmesh circleg Hmax 1 error err 1 while err gt 0 001 p e t refinemesh circleg p e t u assempde circlebl p e t 1 0 1 exact p 1 2 p 2 2 1 4 err norm u exact inf error error err end pdemesh p e t pdesurf p t u pdesurf p t u exact The first MATLAB command creates the initial mesh using the parameterizing function ci rcl eg Examp
110. dary Condition dialog box The Dirichlet condition u x is entered by typingx 2 intothe r edit box Next open the PDE Specification dialog box to define the PDE This is an elliptic equation with ent O Gadi SO 4 1 Vu The nonlinear cis entered into the c edit box as 1 sqrt 1 ux 2 uy 2 Initialize a mesh and refine it once Before solving the PDE select Parameters from the Solve menu and check the Use nonlinear solver option Also set the tolerance parameter to 0 001 Push the button to solve the PDE Use the Plot Selection dialog box to plot the solution in 3 D check u andcontinuous selections in the Height column to visualize the saddle shape of the solution 2 10 Examples of Elliptic Problems Using Command Line Functions Working from the command line the following sequence of commands solves the minimal surface problem and plots the solution The M files ci rcl eg and circleb2 contain the geometry specification and boundary condition functions respectively g circleg b circleb2 gt cal sqrt ltux 2 uy 2 rtol le 3 p e t initmesh g p e t refinemesh g p e t u pdenonlin b p e t c 0 0 Tol Tol rtol pdesurf p t u Y ou can run this example by typingpdedemo3 2 11 2 Examples 2 12 Domain Decomposition The PDE Toolbox is designed to deal with one level domain decomposition If Q has a complicated geometry it is of
111. dge size toi nf inthe Mesh Parameters dialog box Y ou open the dialog box by selecting the Parameters option from the Mesh menu Also select the items Show Node Labels and Show Triangle Labels in the Mesh menu Then create the initial mesh by pressing the A button This can also be done by selecting the Initialize Mesh option from the Mesh menu The figure below appears 0 8F 0 6F 0 4 0 21 inttmesh The corresponding mesh data structures can be exported to the main workspace by selecting the Export Mesh option from the Mesh menu p p 1 1 1 0 0 1 1 1 1 1 0 0 e 1 2 3 4 5 6 2 3 4 5 6 1 0 0 0 0 0 0 1 1 1 1 1 1 1 2 3 4 5 6 1 1 1 1 1 1 0 0 0 0 0 0 t t 1 2 3 1 2 3 4 5 5 5 5 6 1 1 1 1 See Also decsg pdegeom jigglemesh refinemesh Reference George P L Automatic Mesh Generation Application to Finite Element Methods WILEY 1991 jigglemesh Purpose Synopsis Description Examples See Also J iggle internal points of a triangular mesh pl jigglemesh p e t pl jigglemesh p e t PropertyName PropertyValue pl jigglemesh p e t jigglesthetriangular mesh by adjusting the node point positions The quality of the mesh normally increases The following property name property value pairs are allowed Property Value Default Description Opt off mean min mean Optimization method Iter numeric 1 or 20 see below Maximum ite
112. ding edit field is disabled e Plot style contains three pop up menus from which you can control the plot style for the color arrow and height plot types respectively The available plot styles for color surface plots are Interpolated shading A surface plot using the selected colormap and interpolated shading e each triangular area is colored using a linear interpolated shading the default Flat shading A surface plot using the selected colormap and flat shading i e each triangular area is colored using a constant color You can use two different arrow plot styles Proportional The length of the arrow corresponds to the magnitude of the property that you visualize the default Normalized The lengths of all arrows are normalized i e all arrows have the same length This is useful when you are interested in the direction of the vector field The direction is clearly visible even in areas where the magnitude of the field is very small For height 3 D plots the available plot styles are Continuous Produces a smooth continuous plot by interpolating data from triangle midpoints to the mesh nodes the default Discontinuous Produces a discontinuous plot where data and z height are constant on each triangle A total of three properties of the solution two scalar properties and one vector field can be visualized simultaneously If the Height 3 D plot 3 33 3 The Graphical User Interface
113. disjoint sets contains indices of the points that wholly belong to the subdomains listed ins dl andc lists points that also belongs to the other subdomains c pdesdp p e t sdl returns indices of points that belong to more than one of the subdomains in sdi i pdesdt t sdl given triangle datat anda list of subdomain numbers sd i contains indices of the triangles inside that set of subdomains i pdesde e sdl given edge data e it extracts indices of outer boundary edges of the set of subdomains Ifsdl is not given a list of all subdomains is assumed 5 74 pdesmech Purpose Synopsis Description Calculate structural mechanics tensor functions ux pdesmech p t c u PropertyName PropertyValue ux pdesmech p t c u pl vl returns a tensor expression evaluated at the center of each triangle The tensor expressions are stresses and strains for structural mechanics applications with planestress or planestrain conditions pdesmech iSintended to be used for postprocessing of a solution computed using the structural mechanics application modes of the pdet ool GUI after exporting the solution the mesh and the PDE coefficients to the MATLAB main workspace Poisson s ratio nu has to be supplied explicitly for calculations of shear stresses and strains and for the von Mises effective stress in plane strain mode Valid property name property value pairs include Property Name Property Value Default D
114. e AE u R TR 0 dt dt and similarly for H The case of time harmonic fields is treated by using the complex form replacing E by E ca The plane case of this PDE Toolbox mode has Ee 0 0 E and dd 0 0 J 8 and the magnetic field 1 e HH y 0 T xE The scalar equation for E becomes V ive Goc 0 eJE 0 2 51 2 Examples This is the equation used by the PDE Toolbox in the AC Power Electromagnetics mode It is a complex Helmholtz s equation describing the propagation of plane electromagnetic waves in imperfect dielectrics and good conductors o we A complex permittivity can be defined as e jo w The conditions at material interfaces with abrupt changes of e and u arethe natural ones for the variational formulation and need no special attention The PDE parameters that have to be entered into the PDE Specification dialog box are the angular frequency o the magnetic permeability u the conductivity o and the coefficient of dielectricity e The boundary conditions associated with this mode are a Dirichlet boundary condition specifying the value of the electric field E on the boundary and a Neumann condition specifying the normal derivative of E This is equivalent to specifying the tangential component of the magnetic field H Hos Ln QUE Interesting properties that can be computed from the solution the electric field E arethe current density oE and the magnetic flux dens
115. e a growth rate of 30 Hgr ad must be between 1 and 2 Both theBox and nit property are related to the way the mesh algorithm works By turning on Box you can get a good idea of how the mesh generation algorithm works within the bounding box By turning on nit you can seethe initial triangulation of the boundaries By using the command sequence p e t initmesh dl hmax inf init on uxy tn a2 a3 tri2grid p t zeros size p 2 x y n t 4 tn you can determine the subdomain number n of the point x y If the point is outside the geometry tn is NaN and the command n t 4 tn results in a failure The i gg e property is used to control whether jiggling of the mesh should be attempted see the entry onj i ggl emesh for details J iggling can be done until the minimum or the mean of the quality of the triangles decreases i ggl el ter can be used to set an upper limit on the number of iterations 5 39 initmesh Algorithm Examples init mesh implements a Delaunay triangulation algorithm Place node points on the edges Enclose geometry in bounding box Triangulate edges Check that the triangulation respects boundaries Insert node points into centers of circumscribed circles of large triangles Repeat from step 4 if Hmax not yet achieved Remove bounding box YOU A WN P Make a simple triangular mesh of the L shaped membrane in pdet ool Before you do anything in pdet oo set the Maximum e
116. e respectively The rotational angle of the ellipse is stored in row 12 Examples The command sequence below starts pdet oo and draws a unit circleand a unit square pdecirc 0 0 1 pderect 0 1 0 1 Insert the set formula C1 5Q1 Export the Geometry Description matrix set formula and Name Space matrix to the MATLAB main workspace by selecting the Export Geometry Description option from the Draw menu Then type dl bt decsg gd sf ns dl 2 0000 2 0000 1 0000 1 0000 1 0000 0 0 1 0000 0 0000 0 0000 1 0000 0 0 0000 1 0000 1 0000 0 1 0000 0 0000 1 0000 1 0000 0 0 1 0000 0 0 0000 0 0 1 0000 1 0000 1 0000 1 0000 1 0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0000 1 0000 1 0000 bt 1 0 Note that there is one minimal region with five edge segments three circle edge segments and two line edge segments 5 33 decsg Algorithm Cautionary Diagnostics See Also 5 34 The algorithm consists of the following steps 1 Determine the intersection points between the borders of the model objects 2 For each intersection point sort the incoming edge segments on angle and curvature 3 Determine if the induced graph is connected If not add some appropriate edges and redo algorithm from step 1 4 Cycle through edge segments of minimal regions 5 For each original region determine minimal regions inside it 6 Organize output and remove the additional edges The input CSG model is not chec
117. e Export Geometry Description Set Formula Labels option from the Draw menu inpdet ool Each column in the Geometry Description matrix corresponds to an object in the CSG model Four types of solid objects are supported The object type is specified in row 1 e For thecirclesolid row one contains 1 the second and third row contain the center x and y coordinates respectively Row four contains the radius of the circle e For a polygon solid row one contains 2 and the second row contains the number n of line segments in the boundary of the polygon The following n 5 31 decsg 5 32 rows contain the x coordinates of the starting points of the edges and the following n rows contain the y coordinates of the starting points of the edges e For arectanglesolid row one contains 3 The format is otherwise identical to the polygon format e For an dlipsesolid row one contains 4 the second and third row contains the center x and y coordinates respectively Rows four and five contain the semiaxes of the ellipse The rotational angle of the ellipseis stored in row six Set Formula sf contains a set formula expressed with the set of variables listed inns The operators and correspond to the set operations union intersection and set difference respectively The precedence of the operators and are the same has higher precedence The precedence can be controlled with parentheses
118. e PDE Specification dialog box and enter the PDE parameters The angular frequency 2 1 50 y PDE Specification Equation div 1 mu grad E j omega sigma omega 2 epsilon E 0 E electric field Type of PDE Coefficient Value Description 4 Elliptic omega 2 nit50 Angular frequency mu d pit1E 7 Magnetic permeability sigma 57E6 Conductivity epsilon 8 8E 12 Coeff of dielectricity OK Cancel Initialize the mesh and solve the equation Due to the skin effect the current density at the surface of the conductor is much higher than in the conductor s interior This is clearly visualized by plotting the current density J asa 3 D plot To improve the accuracy of the solution close to the surface you need to refine the mesh Open the Solve Parameters dialog box and check the Adaptive mode option Also set the maximum numbers of triangles to nf the maximum numbers of refinements to 1 and use the triangle selection method that picks the worst triangles Recompute the solution several times Each time the adaptive solver refines the area with the largest errors The 2 53 2 Examples number of triangles is printed on the command line The mesh below is the result of successive adaptions and contains 833 triangles z POF Toolbox LUntitled File Edt Options Diaw Boundary PDE Mesh Solve Plot Window Help B mloje gt 39 pe A AA ac Power Electromannetios PE Y 015
119. e Row 2 contains the end parameter value e Row 3 contains the label of the left hand region left with respect to direction induced by start and end from row 1 and 2 e Row 4 contains the label of the right hand region 5 57 pdegeom Examples 5 58 The complement of the union of all regions is assigned the region number 0 x y pdegeom bs s produces coordinates of edge segment points bs specifies the edge segments ands the corresponding parameter values bs can be a scalar The parameter s should be approximately proportional to the curve length All minimal regions should have at least two and preferably three edge segments in their boundary The function cardg defines the geometry of a cardioid r 2 14cos 6 function x y cardg bs s CARDG Geometry File defining the geometry of a cardioid nbs 4 if nargin x nbs return end dl 0 pi 2 pi 3 pi 2 pi 2 pi 3 pi 2 2 pi 1 1 1 1 0 0 0 01 if nargin x dl bs return end x zeros size s y zeros size s m n size bs if m n 1 bs bs ones size s expand bs elseif m size s 1 n size s 2 error bs must be scalar or of same size as s end pdegeom Cautionary See Also nth 400 th linspace 0 2 pi nth r 2 1 cos th xt r cos th yt r sin th th pdearcl th xt yt s 0 2 pi r 2 1 cos th x r cos th y r sin th We use the function pdearc to make the parameter s proportional to arc l
120. e are collected into HU R where H is an M by N matrix and R is an M vector With the reaction force term the system becomes KU H F HU R The constraints can be solved for m of the U variables the remaining called V an Ny m vector The null space of H is spanned by the columns of B and U BV Ug makes U satisfy the Dirichlet conditions A permutation to block diagonal form exploits the sparsity of H to speed up the following computation to find B in a numerically stable way u can be eliminated by premultiplying by B since by the construction HB 0 or B H 0 The reduced system becomes B KBV B F B Kuy which is symmetric and positive definite if K is 4 12 The Parabolic Equation The Parabolic Equation The elliptic solver allows other types of equations to be more or less easily implemented Below we show how the parabolic equation can be reduced to solving elliptic equations This is done by the toolbox function parabolic Consider the equation Y cVu au f in Q with the initial condition u x 0 u x xEQ and boundary conditions of the same kind as for the elliptic equation on dQ The heat equation reads pc V kVu h u u f in the presence of distributed heat loss to the surroundinos p is the density C thermal capacity k thermal conductivity h film coefficient u ambient temperature and f heat source For time independent coefficients the steady state solution of the
121. e tolerances that are passed to the ODE solver ul parabolic u0 tlist K F B ud M produces the solution to the ODE problem du B MB K u F u Bu ug with initial value for u beingu0 Solve the heat equation du oy nH on a square geometry 1 lt x y lt 1 s quar eg Choose u 0 1 on the disk x y lt 0 4 and u 0 0 otherwise Use Dirichlet boundary conditions u 0 squareb1 Compute the solution at times inspace 0 0 1 20 p e t initmesh squareg p e t refinemesh squareg p e t ud zeros size p 2 1 ix find sqrt p 1 2 p 2 2 lt 0 4 ud ix ones size ix tlist linspace 0 0 1 20 ul parabolic u0 tlist squarebl p e t 1 0 1 1 assempde hyperbolic pdeadgsc Purpose Select triangles using a relative tolerance criterion Synopsis bt pdeadgsc p t c a f u errf tol Description bt pdeadgsc p t c a f u errf tol returns indices of triangles to be refinedinbt Used fromadapt mesh to select thetriangles to be further refined The geometry of the PDE problem is given by the mesh data p andt See the entry oni nit mesh for more details c a andf arePDE coefficients Seeassempde for details u is the current solution given as a column vector See the entry on assempde for details errf is the error indicator as calculated by pdej mps tol is a tolerance parameter Triangles are selected using the criterion errf gt tol scale wherescale i
122. e wave equation 2 ou Ay 0 ot for transverse vibrations of a membrane on a square with corners in 1 1 1 1 1 1 and 1 1 The membrane is fixed u 0 at the left and right hand sides and is free o 0 at the upper and lower hand sides Additionally du ty ot boundary conditions for the solution to be well behaved If we start at du O _ ot satisfy the boundary conditions The reason for the arctan and exponential functions is to introduce more modes into the solution ee we need initial values for u ty and Theinitial values need to match the E t 0 u 0 arctan cos 37 and 3 sin x x en are initial values that Using the Graphical User Interface Usethepdetool GUI in the Generic Scalar mode Draw the square usingthe Rectangle square option from the Draw menu or the button with the rectangle icon Proceed to define the boundary conditions by pressing the dQ button and then double click on the boundaries to define the boundary conditions Initialize the mesh by pressing the A button or by selecting Initialize mesh from the Mesh menu 2 23 2 Examples 2 24 Also definethe hyperbolic PDE by openingthe PDE Specification dialog box selecting the hyperbolic PDE and entering the appropriate coefficient values The general hyperbolic PDE is described by 2 d Vv cVu au f ot so for the wave equation you get c 1 a 0 f 0 and d 1 Before solving the PDE select Paramet
123. ecirc xc yc radius label pdecirc xc yc radius draws a circle with center in xc yc and radius radius Ifthepdetool GUI is not active it is automatically started and the circleis drawn in an empty geometry model The optional argument abel assigns a name tothe circle otherwise a default name is chosen The state of the Geometry Description matrix insidepdetoo l is updated to include the circle You can export the Geometry Description matrix from pdetool by using the Export Geometry Description option from the Draw menu The format of the Geometry Description matrix is described in the entry on decsg The command below starts pdet oo and draws a unit cirde pdecirc 0 0 1 pdeellip pdepoly pderect pdetool 5 51 pdecont Purpose Synopsis Description Examples See Also 5 52 Shorthand command for contour plot pdecont p t u pdecont p t u n pdecont p t u v h pdecont p t u h pdecont p t u n h pdecont p t u v pdecont p t u draws 10 level curves of the PDE node or triangle data u h pdecont p t u additionally returns handles to the drawn axes objects Ifu isa column vector node data is assumed If u is a row vector triangle data is assumed Triangle data is converted to node data using the function pdeprtni The geometry of the PDE problem is given by the mesh datap andt Details on the mesh data representation can be found in the entry oni nit mesh pdecont p t u n plots
124. ect to q Row 3 N to 3 N N 1 contain the lengths for the strings representing g Row 3 N N to3 N N MN 1 contain the lengths for the strings representing h The lengths are stored in column wise order with respect to h Row 3 N2 N NM to3 N2 N MN M 1 contain the lengths for the strings representing r The following rows contain MATLAB text expressions representing the actual boundary condition functions The text strings have the lengths according to above The MATLAB text expressions are stored in column wise order with respect to matrices h and q There are no separation characters between the 5 15 assemb Examples 5 16 strings You can insert MATLAB expressions containing the following variables e The 2 D coordinates x and y e A boundary segment parameter s proportional to arc length s is O at the start of the boundary segment and increases to 1 along the boundary segment in the direction indicated by the arrow e Theoutward normal vector components nx andny If you need the tangential vector it can be expressed using nx and ny since tx ny and ty nNx e The solution u only if the input argument u has been specified e Thetimet only if the input argument t i me has been specified The following examples describe the format of the boundary condition matrix For a boundary in a scalar PDE N 1 with Neumann boundary condition M 0 CVu X2 the boundary condition would be re
125. ection method Seep de de mo 7 for an example of a user defined triangle selection method Function parameter The function parameter allows fine tuning of the triangleselection methods F or the worst triangle method pdeadworst it is the fraction of the worst value that is used to determine which triangles to refine For the relative tolerance method it is a tolerance parameter that controls how well the solution fits the PDE Refinement method Can ber egul ar orl ongest Seethe Parameters dialog box description in the Mesh Menu section PDE Toolbox M enus If the problem is nonlinear i e parameters in the PDE are directly dependent on the solution u a nonlinear solver must be used The following parameters are used Use nonlinear solver Toggle the nonlinear solver on off Nonlinear tolerance Tolerance parameter for the nonlinear solver Initial solution An initial guess Can bea constant or a function of x and y given as a MATLAB expression that can be evaluated on the nodes of the current mesh Examples 1 andexp x y Optional parameter defaults to zero J acobian J acobian approximation method f i xed the default a fixed point iteration umped a lumped diagonal approximation or ful the full J acobian Norm Thetypeof norm used for computing the residual Enter asener gy for an energy norm or as a real scalar p to give the p norm The default is Inf the infinity maximu
126. ed geometry by the function i ni t mesh and can be altered by the functions 1 39 1 Tutorial refinemesh andj iggl emesh The Export Mesh option from the Mesh menu provides another way of creating mesh data Theadapt mesh function creates mesh data as part of the solution process The mesh may be plotted with thepdemesh function The details on the mesh data representation are given in the Reference entry oni nit mesh Solution The solution of a PDE problem is represented by the solution vector A solution gives the value at each mesh point of each dependent variable perhaps at several points in time or connected with different eigenvalues Solution vectors are produced fromthe mesh the boundary conditions and the equation coefficients byassempde pdenonlin adaptmesh parabolic hyperbolic and pdeei g The Export Solution option from the Solve menu exports solutions to the workspace Since the meaning of a solution vector is dependent on its corresponding mesh data they are always used together when a solution is presented The details on solution vectors are given in the Chapter 5 Reference entry on assempde Post Processing and Presentation Given a solution mesh pair a variety of tools is provided for the visualization and processing of the data pdei ntrp andpdeprtni can be used to interpolate between functions defined at triangle nodes and functions defined at triangle midpoints tri 2grid interpolates a functions from
127. eee eens 2 27 Using the Graphical User Interface 2 27 Using Command Line Functions 000 ee eee 2 28 L Shaped Membrane with Rounded Corner 2 31 Eigenvalues and Eigenmodes of a Square 0 2 32 Using the Graphical User Interface 2 33 Using Command Line Functions 000 eee 2 33 Application Modes 0 scene eens 2 35 The Application Modes and the GUI 2 35 Structural Mechanics Plane StresS ooooo oooooo 2 36 EXaMpl is ls ee ES 2 39 Using the Graphical User Interface o 2 39 Structural Mechanics Plane Strain o oooooo 2 41 ElectrostaticSiac cit tai soa eee eee aea 2 43 Example score weed wt a theres a tase eae Bede ee SEN 2 44 Using the Graphical User Interface 2 44 Magnetostatics rri aa Soe evista Va cee kee mite ee Wega oes 2 46 Example ici Ree earned AAA 2 47 Using the Graphical User Interface 2 48 AC Power Electromagnetics 00 cece eee eee 2 51 EXA ii aa ape gee ae eke eb 4 2 52 Using the Graphical User Interface 2 53 Conductive Media DG seris troses ian CEEE ARE DEA 2 55 Example 232222440 a Ded e a yee a eS 2 55 Using the Graphical User Interface o o 2 56 Heat Transfer 0 cent e ee enes 2 57 EXAMP nin wind howe Ye E bee ns 2 58 Using the Graphical User Interface
128. eigenvalues in the interval but then many restarts are needed For large values of j max which is the preferred choice mul 1 runs are needed mu is the maximum multiplicity of an eigenvalue in the interval The algorithm works on nonsymmetric as well as symmetric pencils but then accuracy is approximately to times the Henrici departure from normality The parameter s pd is used only to choose between s y mmmd and col mmd when factorizing the former being marginally better for symmetric matrices close to the lower end of the spectrum In case of trouble e f convergence is too slow try in this order of priority 1 asmaller interval b ub 2 alargerjmax 3 alarger max mul e f factorization fails try again with b or ub finite Then shift is chosen at random and hopefully not at an eigenvalue If it fails again check whether pencil may be singular e fit goes on forever there may be too many eigenvalues in thestrip Try with a small value ma x mul 2 and see which eigenvalues you get Those you get are 5 93 sptarn See Also Reference 5 94 some of the eigenvalues but a negativei result tells you that you have not gotten them all e f memory overflow try smaller j max e The algorithm is designed for eigenvalues close to the real axis If you want those close to the imaginary axis try A i A e When spd the shift is at b sothat advantage is taken of the faster factorization for symmetric positive def
129. ength You can test the function by typing pdegplot cardg axis equal p e t initmesh cardg pdemesh p e t axis equal Then solve the PDE problem Au 1 on the geometry defined by the the cardioid Use Dirichlet boundary conditions u 0 on dQ Finally plot the solution U assempde cardb p e t 1 0 1 pdesurf p t u The parameter s should be approximately proportional to the curve length All minimal regions should have at least two and preferably three edge segments in their boundary initmesh refinemesh pdearcl 5 59 pdegplot Purpose Synopsis Description Examples See Also 5 60 Plot a PDE geometry pdegplot g h pdegplot g pdegplot g plots the geometry of a PDE problem h pdegpl ot 9 returns handles to the plotted axes objects y describes the geometry of the PDE problem y can either be a Decomposed Geometry matrix or the name of a Geometry M file The formats of the Decomposed Geometry matrix and Geometry M file are described in the entries ondecsg and pdegeom respectively Plot the geometry of the L shaped membrane pdegplot shapeg pdegeom pdegrad Purpose Synopsis Description See Also The gradient of a PDE solution ux uy pdegrad p t u ux uy pdegrad p t u sdl ux uy pdegrad p t u returns the gradient of u evaluated at the center of each triangle Row i from 1toN of ux contains du Ox row i from 1toN of uy co
130. enmodes decrease as the radius increases and the shape of the L shaped membrane becomes more rounded In the figure below the first eigenmode of an L shaped membrane with a rounded corner is plotted 2 31 2 Examples 2 32 vf PDE Toolbox LSHAPEC M File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help X 1 132 Y 0 9151 Info Specify type of PDE solution plot Exit m The first eigenmode for an L shaped membranewith a rounded corner Eigenvalues and Eigenmodes of a Square Let s study the eigenvalues and eigenmodes of a square with an interesting set of boundary conditions The square has corners in 1 1 1 1 1 1 and 1 1 The boundary conditions are as follows e On the left boundary the Dirichlet condition u 0 e On the upper and lower boundary the Neumann condition eu 0 e On theright boundary the generalized Neumann condition ee 2u 0 n The eigenvalue PDE problem is Au u We are interested in the eigenvalues smaller than 10 and the corresponding eigenmodes so the search rangeis Inf 10 Thesign in the generalized Neumann condition is such that there are negative eigenvalues Examples of Eigenvalue Problems Using the Graphical User Interface Usingthepdetool GUI in the Generic Scalar mode draw the square using the Rectangle square option from the Draw menu or the button with the rectangle icon
131. entation oo ooooo 1 40 Hints and Suggestions for Using Command Line FUNCUONS ars tite dad Olt da dl oak 1 40 Examples 2 Examples of Elliptic Problems 2 2 Poisson s Equation on Unit Disk 0 00 cece eee 2 2 Using the Graphical User Interface oooo 2 2 Using Command Line Functions 0 eee eee 2 4 A ScatteringProbleM o ooocoocoocc ee 2 6 Using the Graphical User Interface o o o o 2 8 A Minimal Surface Problem 000 c eee eee eens 2 10 Using the Graphical User Interface 2 10 Using Command Line Functions 000 eee ee 2 11 Domain Decomposition 0 0 0 cece eee 2 12 Examples of Parabolic ProblemsS 2 16 The Heat Equation A Heated Metal Block 2 16 Using the Graphical User Interface 2 17 Using Command Line Functions 000 eee ee 2 19 Heat Distribution in Radioactive Rod 2 21 Using the Graphical User Interface 2 22 Examples of Hyperbolic Problems 2 23 The Wave Equation 00 0 eee ee 2 23 Using the Graphical User Interface 2 23 Using Command Line Functions 00 ee eee 2 25 ii Contents Examples of Eigenvalue Problems 2 27 Eigenvalues and Eigenfunctions for the L Shaped Membrane 0 0 0 ee eee
132. equation is the solution to our standard elliptic equation V cvu au f Assuming a triangular mesh on Q and at any timet 0 expand the solution to the PDE as a function of x in the Finite Element method basis u x t Y U 00 00 4 13 4 The Finite Element M ethod Plugging the expansion into the PDE multiplying with a test function 9 integrating over Q and applying Green s formula and the boundary conditions yield dU t De at D a cVo ad dx 29 9 ds U LO Sa 08 vj In matrix notation we have to solve the linear large and sparse ODE system MC KU dt This method is traditionally called method of lines semidiscretization Solving the ODE with the initial value U 0 u x yields the solution tothe PDE at each node x and timet Note that K andF are the stiffness matrix and the right hand side of the elliptic problem V Vu au f inQ with the original boundary conditions while M is just the mass matrix of the problem V 0Vu du 0 in Q When the Dirichlet conditions are time dependent F contains contributions from time derivatives of h and r These derivatives are evaluated by finite differences of the user specified data 4 14 The Parabolic Equation The ODE system is ill conditioned Explicit time integrators are forced by stability requirements to very short time steps while implicit solvers can be expensive since they solve an elliptic problem at ev
133. ercial computer software or commercial software documentation was obtained as set forth in subparagraph a of the Rights in Commercial Computer Software or Commercial Software Documentation Clause at DFARS 227 7202 3 therefore the rights set forth herein shall apply and b For any other unit or agency NOTICE Notwithstanding any other lease or license agreement that may pertain to or accompany the delivery of the computer software and accompanying documentation the rights of the Government regarding its use reproduction and disdo sure are as set forth in Clause 52 227 19 c 2 of the FAR MATLAB Simulink Handle Graphics and Real Time Workshop are registered trademarks and Stateflow and Target Language Compiler are trademarks of The MathWorks Inc Other product or brand names are trademarks or registered trademarks of their respective holders Printing History August 1995 First printing February 1996 Reprint Contents Tutorial 1 Introduction oooocco eens 1 2 What Does this Toolbox Do cee ee 1 2 Can UsethePDE Toolbox naana cee ee 1 2 What Problems Can Solve 0 cece ees 1 3 In Which Areas Can the Toolbox Be Used 1 5 How Dol DefineaPDE Problem 0 eae eee 1 5 How Can Solve a PDE Problem 000 cence eee 1 6 Can Use the Toolbox for Nonstandard Problems 1 6 How Can Visualize My Results 02 eee eee 1 6 Are
134. erged or else with the best unconverged approximate eigenvector yj In each step j of this second Arnoldi run the vector is made orthogonal to all vectors in V including the converged Schur vectors from the previous runs This way the algorithm is applied to a projected matrix and picks up a second copy of any double eigenvalue there may be in the interval If anything in the interval converges during this second run athird is attempted and so on until no more approximate eigenvalues 6 show up inside Then the algorithm signals convergence If there are still unconverged approximate eigenvalues after a prescribed maximum number of steps the algorithm signals nonconvergence and reports all solutions it has found 4 19 4 The Finite Element M ethod This is a heuristic strategy that has worked well on both symmetric nonsymmetric and even defective eigenvalue problems There is a tiny theoretical chance of missing an eigenvalue if all the random starting vectors happen to be orthogonal to its eigenvector Normally the algorithm restarts p times if the maximum multiplicity of an eigenvalueis p At each restart a new random starting direction is introduced Note that the shifted and inverted matrix A K uM M is needed only to operate on a vector vj in the Arnoldi algorithm This is done by computing an LU factorization P K uM Q LU using the sparse MATLAB command u P and Q are permutations that make thetriangular factor
135. ers from the Solve menu to open the Solve Parameters dialog box As a list of times enter i nspace 0 5 31 and as initial values for u atan cos pi 2 x ou and for gt 3 sin pi x exp sin pi 2 y EA Solve Parameters Time linspace 0 5 31 u tO atan cos pif2 x uito 3 sinpitid exp sin pir2 y Relative tolerance 0 01 Absolute tolerance 0 001 OK Cancel Examples of Hyperbolic Problems Finally press the button to compute the solution The best plot for viewing the waves moving in the x and y directions is an animation of the whole sequence of solutions Animation is a very real time and memory consuming feature so you may have to cut down on the number of times at which to compute a solution A good suggestion is to check the Plot in x y grid option Using an x y grid can speed up the animation process significantly Using Command Line Functions From the command line solve the equation with the boundary conditions and the initial values above starting at time O and then every 0 05 seconds for five seconds The geometry is described in thefiles quar eg m and the boundary conditions in the files quareb3 m The following sequence of commands then generates a solution and animates it First create a mesh and define the initial values and the times for which you want to solve the equation p e t initmesh squareg SED Tee ye ap ice 4
136. ery time step The numerical integration of the ODE system is performed by the MATLAB ODE Suite functions which are efficient for this class of problems The time step is controlled to satisfy a tolerance on the error and factorizations of coefficient matrices are performed only when necessary When coefficients are time dependent the necessity of re evaluating and re factorizing the matrices each time step may still make the solution time consuming although parabolic re evaluates only that which varies with time In certain cases a time dependent Dirichlet matrix h t may cause the error control to fail even if the problem is mathematically sound and the solution u t is smooth This can happen because the ODE integrator looks only at the reduced solution v with u Bv ud As h changes the pivoting scheme employed for numerical stability may change the elimination order from one step to the next This means that B v and ud all change discontinuously although u itself does not 4 15 4 The Finite Element M ethod The Hyperbolic Equation Using the same ideas as for the parabolic equation hyperbolic implements the numerical solution of 2 dou V cVu tau f in Q ot with the initial conditions u x 0 uy x Mes 0 vox xEQ and usual boundary conditions In particular solutions of the equation Ug CAU 0 are waves moving with speed c Using a given triangulation of Q the method of lines yields the second order ODE s
137. escription tensor ux uy vx vy exx eyy exy sxx syy sxy el Tensor expression e2 s1 52 Yon Mises application ps Hopn Plane stress plane strain nu scalar or string expression 3 Poisson s ratio The available tensor expressions are es yy Buyer a y gt e exx the x direction strain e eyy the y direction strain ey exy the shear strain yyy 5 75 pdesmech Examples 5 76 sxx the x direction stress ox syy the y direction stress oy sxy the shear stress tyy el the first principal strain e4 e2 the second principal strain e2 51 the first principal stress 01 2 the second principal stress o gt 2 2 von Mises the von Mises effective stress o 0 0 0 for plane stress conditions or 0 1 0 0 2v 2v 1 for plane strain conditions Assuming that a problem has been solved using the application mode Structural Mechanics Plane Stress discussed in Structural Mechanics Plane Stress on page 2 36 and that the solution u the mesh datap andt and the PDE coefficient c all have been exported to the MATLAB main workspace the x direction strain is computed as sx pdesmech p t c u tensor sxx To compute the von Mises effective stress for a plane strain problem with Poisson s ratio equal to 0 3 type mises pdesmech p t c u tensor von Mises application pn nu 0 3 pdesurf Purpose Synopsis Description E
138. esponding eigenfunctions p e t initmesh squareg p e t refinemesh squareg p e t The eigenvalue PDE coefficients c a and d for this problem are c 1 a 0 and d 1 You can enter the eigenvalue ranger as the MATLAB vector Inf 10 pdeei g returns two output arguments the eigenvalues as an array anda matrix v of corresponding eigenfunctions v 1 pdeeig squareb2 p e t 1 0 1 Inf 10 2 33 2 Examples To plot the fourth eigenfunction as a surface plot type gt pdesurf p t v 4 This problem is separable i e u x y f x g y The functions f and g are eigenfunctions in the x and y directions respectively In the x direction the first eigenmode is a slowly increasing exponential function The higher modes include sinusoids In the y direction the first eigenmode is a straight line constant the second is half a cosine the third is a full cosine the fourth is one and a half full cosines etc These eigenmodes in the y direction are associated with the eigenvalues A 4 4 There are five eigenvalues smaller than 10 for this problem and the first one is even negative 0 4145 It is possible to trace the eigenvalues above in the eigenvalues of the solution Looking at a plot of the first eigenmode you can see that it is made up of the first eigenmodes in the x and y directions The second eigenmode is made up of the first eigenmode in the x direction and the second eigenmode in the
139. etic permeability Since V B 0 there exists a magnetic vector potential A such that B VxA and Vx iy xA J u The plane case assumes that the current flows are parallel tothez axis so only the z component of A is present A 0 0 A J 0 0 and the equation above can be simplified to the scalar elliptic PDE Y Eva J u where J x y For the 2 D case we can compute the magnetic flux density B as _ aA A oo 55 300 and the magnetic field H in turn is given by H 5B u 2 46 Application Modes The interface condition across subdomain borders between regions of different material properties is that H x n be continuous This implies the continuity of I and does not require special treatment since we are using the variational formulation of the PDE problem In ferromagnetic materials u is usually dependent on the field strength B VA so the nonlinear solver is needed The Dirichlet boundary condition specifies the value of the magnetostatic potential A on the boundary The Neumann condition specifies the value of the normal component of n IVA on the boundary This is equivalent to specifying the tangential value of the magnetic field H on the boundary Visualization of the magnetostatic potential A the magnetic field H and the magnetic flux density B is available B and H can be plotted as vector fields See 2 for a more detailed discussion on Maxwell s equations and magn
140. etostatics Example Asan example of a problem in magnetostatics consider determining the static magnetic field due to the stator windings in a two pole electric motor The motor is considered to be long and when end effects are neglected a 2 D computational model suffices The domain consists of four regions e Two ferromagnetic pieces the stator and the rotor e The air gap between the stator and the rotor e The armature coil carrying the DC current The magnetic permeability u is Lin the air and in the coil In thestator and the rotor u is defined by Umax 2 1 cellV A min Umax 5000 umin 200 and c 0 05 are values that could represent transformer steel The current density J is O everywhere except in the coil where it is 1 2 47 2 Examples 2 48 The geometry of the problem makes the magnetic vector potential A symmetric with respect to y and anti symmetric with respect to x so you can limit the domain to x 0 y20 with the Neumann boundary condition n VA Oonthe x axis and the Dirichlet boundary condition A 0 on the y axis The field outside the motor is neglected leading to the Dirichlet boundary condition A 0 on the exterior boundary Using the Graphical User Interface The geometry is complex involving five circular arcs and two rectangles Using thepdetool GUI set the x axis limits to 1 5 1 5 and the y axis limits to 1 1 Set the application mode to Magnetostatics and use a grid spac
141. ewton iteration tends to the minimizer of the residual i e the solution of miny p U 4 22 N onlinear Equations Closely related to the above problem is the choice of the initial guess U By default the solver sets U and then assembles the FEM matrices K and F and computes UY KV7F The damped Gauss N ewton iteration is then started with U which should be a better guess than U If the boundary conditions do not depend on the solution u then U satisfies them even if U does not Furthermore if the equation is linear then U is the exact FEM solution and the solver does not enter the Gauss N ewton loop There are situations where UV 0 makes no sense or convergence is impossible In some situations you may already have a good approximation and the nonlinear solver can be started with it avoiding the slow convergence regime This idea is used in the adaptive mesh generator It computes a solution U on a mesh evaluates the error and may refine certain triangles The interpolant of U is avery good starting guess for the solution on the refined mesh In general the exact J acobian y 20 n JU is not available Approximation of J n by finite differences in the following way is expensive but feasible The ith column of J n can be approximated by AUD 29 p U which implies the assembling of the FEM matrices for thetriangles containing grid point i A very simple approximation to p which gives a
142. gative gradient of u c Vu user entry A MATLAB expression px py returning a 2 by ntri matrix of data defined on the triangles of the current triangular mesh ntri is the number of triangles in the current mesh The solution u its derivatives u x and uy the x and y components of c Vu cux andcuy andx andy areall availablein the local workspace Data defined on the nodes is interpolated to triangle centers Y ou enter the expression into the edit field tothe right of the Property pop up menu in the User entry column Examples ux uy x y PDE Toolbox M enus For the generic system case the properties available for visualization using color contour lines or z height are u v abs u v and a user entry For visualization using arrows or a deformed mesh you can choose u v or a user entry For applications in structural mechanics u and v are the x and y displacements respectively For the visualization options in the other application modes see Application Modes on page 2 35 Note that the variables available in the local workspace for a user entered expression are the same for all scalar and system modes the solution is always referred to as u and in the system case v e User entry contains four edit fields where you can enter your own expression if you select the user entry property from the corresponding pop up menu to the left of the edit fields If the user entry property is not selected the correspon
143. ge and open a dialog box by double clicking on the boundary or by using the Specify Boundary Conditions option from the Boundary menu 4 Initializethetriangular mesh Press the A button or usethe corresponding Mesh menu option Initialize Mesh Normally the mesh algorithm s default parameters generate a good mesh If necessary they can be accessed using the Parameters menu item 5 If you need a finer mesh the mesh can be refined by pressing the Refine button Pressing the button several times causes a successive refinement of the mesh The cost of a very fine mesh is a significant increase in the number of points where the PDE is solved and consequently a significant increase in the time required to compute the solution Don t refine unless it is required to achieve the desired accuracy For each refinement the number of triangles increases by a factor of four A better way to increase the 1 32 Using the Graphical User Interface accuracy of the solution to elliptic PDE problems is to use the adaptive solver which refines the mesh in the areas where the estimated error of the solution is largest See the entry on adapt mesh in Chapter 5 Reference for an example of how the adaptive solver can solve a Laplace equation with an accuracy that requires more than 10 times as many triangles when regular refinement is used 6 Specify the PDE from the PDE Specification dialog box You can access that dialog box using the PDE bu
144. ge is joined with the the opposing corner and with the other midpoint If only the longest edge is divided its midpoint is joined with the opposing corner initmesh pdegeom pdesdt pdeent 5 91 sptarn Purpose Synopsis Description 5 92 Solve generalized sparse eigenvalue problem xv lmb iresult sptarn A B 1b ub xv mb iresult sptarn A B b ub spd xv mb iresult sptarn A B 1b ub spd tolconv xv mb iresult sptarn A B lb ub spd tolconv j max xv mb iresult sptarn A B b ub spd tolconv j max maxmul xv mb iresult sptarn A B lb ub spd tolconv j max maxmul finds eigenvalues of the penai A AB x 0 in interval Ib ub A matrix of linear polynomials Aj ABj A AB is called a pencil A andB are sparse matrices b and ub are lower and upper bounds for eigenvalues to be sought We may have b i nf if all eigenvalues to the left of ub are sought andr b i nf if all eigenvalues to the right of b are sought One of b andub must be finite A narrower interval makes the algorithm faster In the complex case the real parts of Imb are compared to b and ub xv are eigenvectors ordered sothat nor m a xv b xv diag mb issmall mb is the sorted eigenvalues Ifi resul t gt 0 the algorithm succeeded and all eigenvalues in the intervals have been found Ifi resul t lt 0 the algorithm has not yet been successful there may be more eigenvalues
145. ge to the mesh by selecting the Mesh menu option Undo Mesh Change Initialize the mesh then refine it once and finally jiggle it once x PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help A OS gt Rr AJA AM A Generic Saar HELE ER Set formula Hie ALCA 1 F T T T T 0 8 H KINNE SL 06 p Vek KES SNE AAA ROK se OS 7 X Pe 1 Y x iS S PAA TA VAVA VATAN NM x Kt VAY AVAVAN VAV va Ki Es RE KR WA D 1 15 1 0 5 a as 1 15 Info Click and drag at center to create ellipse Exit 1 Wearenow ready tosolvethe problem Press the button or select Solve PDE from the Solve menu tosolvetheP DE Thesolution is then plotted By default the plot uses interpolated coloring and a linear color map A colorbar is also provided to map the different shades tothe numerical values of the solution If you want the solution can be exported as a vector to the MATLAB main workspace 1 15 1 Tutorial 1 16 x__ _ _ _ ____ RE v PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help i i i 15 Info Push to refine mesh Exit 1 There are many more plot modes available to help you visualize the solution Press the button with the 3 D solution
146. ght part of the dialog box contains a pop up menu with all eigenvalues The plotted solution is the eigenvector associated with the selected eigenvalue By default the smallest eigenvalue is selected You can rotate the 3 D plots by clicking on the plot and while keeping the mouse button down moving the mouse For guidance a surrounding box appears When releasing the mouse the plot is redrawn using the new viewpoint Initially the solution is plotted using 37 5 degrees horizontal rotation and 30 degrees elevation If you press the Plot button the solution is plotted immediately using the current plot setup If there is no current solution available the PDE is first solved The new solution is then plotted The dialog box remains on the screen If you press the Done button the dialog box is closed The current setup is saved but no additional plotting takes place If you press the Cancel button the dialog box is closed The setup remains unchanged since the last plot PDE Toolbox M enus Window Menu From the Window menu you can select all currently open MATLAB figure windows The selected window is brought to the front Help Menu Help Display a brief help window About Display a window with some program information 3 37 3 The Graphical User Interface 3 38 The Toolbar The toolbar underneath the main menu at the top of the GUI contains icon buttons that provide quick and easy access to some of the mos
147. h description tothe MATLAB workspace From the command line or M files you can call functions from the toolbox to do the hard work e g generate meshes discretize your problem perform interpolation plot data on unstructured grids etc while you retain full control over the global numerical algorithm Intro duction What Problems Can I Solve The basic equation of the PDE Toolbox is the PDE V cVu au f ing which we shall refer to as the dliptic equation regardless of whether its coefficients and boundary conditions make the PDE problem elliptic in the mathematical sense Analogously we shall use the terms parabolic equation and hyperbolic equation for equations with spatial operators like the one above and first and second order time derivatives respectively Q is a bounded domain in the plane c a f and the unknown u are scalar complex valued functions defined on Q ccan bea 2 by 2 matrix function on Q The toolbox can also handle the parabolic PDE Ly cVu au f the hyperbolic PDE 3 d v cVu au f 2 ot and the eigenvalue problem V cVu au Adu where d is a complex valued function on Q and is an unknown eigenvalue For the parabolic and hyperbolic PDE the coefficients c a f and d can depend on time A nonlinear solver is available for the nonlinear elliptic PDE V c u Vu a u u fu wherec a and f arefunctions ofthe unknown solution u All solvers can handle the system case
148. hand side Neumann condition n bu 0on all other boundaries Neumann condition n Also for the heat equation we need an initial value the temperature in the metal block at the starting time t In this case the temperature of the block is 0 degrees at the time we start applying heat Finally to complete the problem formulation we specify that the starting time is O and that we want to study the heat distribution during the first five seconds 2 16 Examples of Parabolic Problems Using the Graphical User Interface Once you have started the pdet oo GUI and selected the Generic Scalar mode drawing the CSG model can be done very quickly Draw a rectangle R1 with thecornersinx 0 5 0 5 0 5 0 5 andy 0 8 0 8 0 8 0 8 Draw another rectangle R2 to represent the rectangular cavity Its corners should have the coordinates x 0 05 0 05 0 05 0 05 and y 0 4 0 4 0 4 0 4 Toassist in drawing the narrow rectangle representing the cavity open the Grid Spacing dialog box from the Options and enter x axis extra ticks at 0 05 and0 05 Then turn on the grid and the snap to grid feature A rectangular cavity with the correct dimensions is then easy to draw eee E Grid Spacing x axis linear spacing Auto 1 5 0 5 1 5 x axis extra ticks 0 05 0 05 Y axis linear spacing E Auto Y axis extra ticks Apply Done The CSG model of the metal block is now simply expres
149. hange by clicking to select one boundary segment by Shift clicking to select multi ple segments or by using the Edit menu option Select All to select all boundary segments The selected boundary segments are indicated by black color 1 12 G etting Started For this problem change the boundary condition for all the circle arcs Select them by using the mouse and Shift click on those boundary segments wifi PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help X 1 5 eal Info Click and drag at center to create rectangle Double clicking anywhere on the selected boundary segments opens the Boundary Condition dialog box Here you select the type of boundary condition and enter the boundary condition as a MATLAB expression Change the boundary condition along the selected boundaries to a Neumann condition dn du 5 This means that the solution has a slope of 5 in the normal direction for these boundary segments In the Boundary Condition dialog box select the Neumann condition type and enter 5 in the edit box for the boundary condition parameter g To define a pure Neumann condition leave theq parameter at its default value 0 When you press the OK button notice how the selected boundary segments change to blue to indicate Neumann boundary condition 1 13 1 Tutorial J Next specify the P
150. he boundary conditions 3 Definethe PDE 4 Create the triangular mesh 5 Solvethe PDE 6 Plot the solution and other physical properties calculated from the solution post processing Thepdet ool GUI is designed in asimilar way Y ou work in six different modes each corresponding to one of the steps in the PDE solving process e In Draw mode you can create the 2 D geometry using the constructive solid geometry CSG model paradigm A set of solid objects rectangle circle ellipse and polygon is provided These objects can be combined using set formulas in a flexible way e n Boundary mode you can specify the boundary conditions You can have different types of boundary conditions on different boundaries n this mode the original shapes of the solid objects constitute borders between subdomains of the model Such borders can be eliminated in this mode 1 26 Using the Graphical User Interface e n PDE mode you can interactively specify the type of PDE problem andthe PDE coefficients You can specify the coefficients for each subdomain independently This makes it easy to specify e g various material properties in a PDE model In Mesh mode you can control the automated mesh generation and plot the mesh In Solve mode you can invoke and control the nonlinear and adaptive solver for elliptic problems F or parabolic and hyperbolic PDE problems you can specify the initial values and the times for which the output
151. hese properties can be selected from pop up menus in the Plot Selection dialog box A combination of scalar and vector properties can be plotted simultaneously by selecting different properties to be represented by color height vector field arrows and displacements in a 3 D plot 2 40 Application M odes Select to plot the von Mises effective stress using color and the displacement vector field u v using a deformed mesh Check the Color and Deformed mesh plot types To plot the von Mises effective stress select von Mises from the pop up menu in the Color row In areas where the gradient of the solution the stress is large you need to refine the mesh in order to increase the accuracy of the solution Select Parameters from the Solve menu and check the Adaptive mode option You can use the default options for adaption which are the Worst triangles triangle selection method with the Worst triangle fraction set to 0 5 Now solve the plane stress problem again Check the Show Mesh option in the Plot Selection dialog box to see how the mesh is refined in areas where the stress is large PDE Toolbox PLATE M zi File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help ii D gt r Aa Set formula X 0 3333 Y 0 3333 A R structural Mech Plane Stress Color von Mises Displacement uv T T L 1 L 0 0 3333 0 6667 1 Into Specify type of PD
152. indicates that the ellipse has positive semiaxes 1 indicates that at least one of the semiaxes is nonpositive and 2 indicates that the ellipse is not unique Ifgstat consists of zero entries only then gd is valid and can be used as input argument by decsg decsg csgdel Purpose Synopsis Description See Also Delete borders between minimal regions dl 1 bt1 csgdel dl bt bl dl 1 bt1 csgdel dl bt dl 1 bt1 csgdel dl bt bl deletes the border segmentsinthelistb Ifthe consistency of the Decomposed Geometry matrix is not preserved by deleting the elements in the list b additional border segments are deleted Boundary segments cannot be deleted Theconcepts or border segments boundary segments and minimal regions are explained in the entry decsg dl and di 1 are Decomposed Geometry matrices See the entry ondecsg fora description of the Decomposed Geometry matrix The format of the Boolean tables bt and bt 1 is also described in the entry ondecsg dl 1 bt1 csgdel dl bt deletes all border segments decsg csgchk 5 29 decsg Purpose Synopsis Description 5 30 Decompose Constructive Solid Geometry into minimal regions dl decsg gd dl decsg gd sf ns dl bt decsg gd dl bt decsg gd sf ns dl bt dl1 bt1 msb decsg gd dl bt dl1 bt1 msb decsg gd sf ns This function analyzes the Constructive Solid Geometry modd CSG model that you draw It an
153. ing of 0 1 The model is a union of circles and rectangles the reduction tothe first quadrant is achieved by intersection with a square Using the snap to grid feature you can draw the geometry using the mouse or you can draw it by entering the following commands pdecirc 0 0 1 C1 pdecirc 0 0 0 8 C2 pdecirc 0 0 0 6 C3 pdecirc 0 0 0 5 C4 pdecirc 0 0 0 4 C5 pderect 0 2 0 2 0 2 0 9 R1 pderect 0 1 0 1 0 2 0 9 R2 pderect 0 1 0 1 50Q1 Y ou should get a CSG model similar to the one in the plot below fa PDE Toolbox MAGNET M File Edil Options Draw Boundary PDE Mesh Solve Plot Window Help Info Click and drag at center to create rectangle Exit Application M odes Enter the following set formula to reduce the model to the first quadrant C1 C2 C3 C4 C5 R1 R2 SQ1 In boundary mode you need to remove a number of subdomain borders Using Shift click select borders and remove them using the Remove Subdomain Border option from the Boundary menu until the geometry consists of four subdomains the stator the coil the air gap and the rotor In the plot below the stator is subdomain 1 the rotor is subdomain 2 the coil is subdomain 3 and the air gap is subdomain 4 Note that the numbering of your subdomains may be different Az PDE Toolbox MAGNET M File Edit Options Draw Boundary PDE Mesh So
154. ing boundary condition types On a generalized Neumann boundary segment q and g are related tothe normal derivative value by cVu qu g On a Dirichlet boundary segment hu r The toolbox can also handle systems of partial differential equations over the domain Q Let the number of variables in the system be N Our general boundary condition is hu r Vu qu g thy By the notation c Vu we mean the N by 1 matrix with i 1 component assemb N d t on iss d coste kiaz cos 0 c 1 27 sin O c dine sin Q c j 2 ih where a is the angle of the normal vector of the boundary pointing in the direction out from Q the domain The Boundary Condition matrix is created internally in pdet oo actually a function called by pdet 00 and then used from the function as semb for assembling the contributions from the boundary to the matrices Q G H and R The Boundary Condition matrix can also be saved onto a file as a boundary M file for later use with the wbound function For each column in the Decomposed Geometry matrix there must bea corresponding column in the Boundary Condition matrix The format of each column is according to the list below Row one contains the dimension N of thesystem Row two contains the number M of Dirichlet boundary conditions Row three to 3 N 1 contain the lengths for the strings representing q The lengths are stored in column wise order with resp
155. inite matrices No harm is done but the execution is slower if b is above the lowest eigenvalue pdeei g Golub Gene H and Charles F Van Loan Matrix Computations 2nd edition Johns Hopkins University Press Baltimore MD 1989 Saad Y ousef Variations on Arnoldi s Method for Computing Eigenelements of Large Unsymmetric Matrices Linear Algebra and its Applications Vol 34 1980 pp 269 295 tri2grid Purpose Synopsis Description See Also Interpolate from PDE triangular mesh to rectangular grid uxy tri2grid p t u x y uxy tn a2 a3 tri2grid p t u x y uxy tri2grid p t u tn a2 a3 uxy tri2grid p t u x y computes the function values uxy over the grid defined by the vectorsx and y fromthefunction u with values on thetriangular mesh defined byp and t Values are computed using linear interpolation in the triangle containing the grid point The vectors x and y must be increasing uxy tn a2 a3 tri2grid p t u x y additionally lists the indextn of the triangle containing each grid point and interpolation coefficients a2 and a3 uxy tri2grid p t u tn a2 a3 Wwithtn a2 anda3 computed in an earlier call totri2grid interpolates using the same grid as in the earlier call This variant is however much faster if several functions have to be interpolated using the same grid For grid points outside of the triangular mesh NaN is returned inuxy tn a2 and a3 initmesh refinemesh assempde 5 95 w b
156. iption Examples See Also Draw rectangle pderect xy pderect xy label pderect xy draws a rectangle with corner coordinates defined by x y x mi n xmax ymin ymax Ifthepdetool GUI is not active it is automatically started and the rectangle is drawn in an empty geometry model The optional argument abel assigns a name to the rectangle otherwise a default name is chosen The state of the Geometry Description matrix inside pdet oo is updated to include the rectangle You can export the Geometry Description matrix from pdetool by selecting the Export Geometry Description option from the Draw menu Theformat of the Geometry Description matrix is described in the entry on decsg The command sequence below starts pdetool and draws the L shaped membrane as the union of three squares pderect 1 0 1 0 pderect 0 1 1 0 pderect 0 1 0 1 pdecirc pdeellip pdepoly pdetool 5 73 pdesdp pdesde pdesdt Purpose Indices of points edges triangles in a set of subdomains Synopsis c pdesdp p e t i cl pdesdp p e t c pdesdp p e t sdl i C pdesdp p e t sdl i pdesdt t i pdesdt t sdl i pdesde e i pdesde e sdl Description i c pdesdp p e t sdl given mesh datap e andt anda list of subdomain numbers sd the function returns all points belonging to those subdomains A point can belong to several subdomains and the points belonging to the domains insd are divided into two
157. iption Triangle geometry data ar al a2 a3 pdetrg p t ar glx gly g2x g2y 93x g3y pdetrg p t ar al a2 a3 pdetrg p t returns thearea of each triangleinar and half of the negative cotangent of each angle inal a2 anda3 ar glx gly g2x g2y 93x g3y pdetrg p t returns the area and the gradient components of the triangle base functions The triangular mesh of the PDE problem is given by the mesh datap andt Details on the mesh data representation can be found in the entry oni ni t mesh 5 83 pdetriq Purpose Synopsis Description See Also Reference 5 84 Triangle quality measure q pdetriq p t q pdetriq p t returns a triangle quality measure given mesh data Thetriangular mesh is given by the mesh data p e andt Details on the mesh data representation can be found in the entry on i nit mesh The triangle quality is given by the formula _ 48a TER ae where a is the area and hy hz and h3 the side lengths of the triangle If q gt 0 6 the triangle is of acceptable quality q 1 when h h h3 initmesh jigglemesh refinemesh Bank Randolph E PLTMG A Software Packagefor Solving Elliptic Partial Differential Equations User s Guide 6 0 Society for Industrial and Applied Mathematics Philadelphia PA 1990 poiasma Purpose Synopsis Description See Also Boundary point matrix contributions for fast solvers of Poisson s equation K poi asma nl n2 hl h2 K poiasma
158. it to 0 In the diamond shaped region enter a density of 1 a heat capacity of 0 1 and a coefficient of heat conduction of 2 Enter 4 in the edit field for the heat source Thetransversal heat transfer term h Tax T is not used so set h the convective heat transfer coefficient to 0 Since you are solving a dynamic PDE you have to define an initial value and the times at which you want to solvethe PDE Open the Solve Parameters dialog box by selecting Parameters from the Solve menu The dynamics for this problem is very fast the temperature reaches steady state in about 0 1 time units To capture the interesting part of the dynamics enter ogs pace 2 1 10 as the vector of times at which to solve the heat equation logspace 2 1 10 gives 10 logarithmically spaced numbers between 0 01 and 0 1 Set theinitial value of the temperature to 0 If the boundary conditions and the initial value differ the problem formulation contains discontinuities Solve the PDE By default the temperature distribution at the last time is plotted The best way to visualize the dynamic behavior of the temperature is to animate the solution When animating turn on the Height 3 D plot option to animate a 3 D plot Also select the Plot in x y grid option Usinga rectangular grid instead of a triangular grid speeds up the animation process significantly 2 59 2 Examples Other interesting visualizations are made by plotting isothermal lines
159. ity B LVxE 0 The electric field E the current density J the magnetic field H and the magnetic flux density B are available for plots Additionally the resistive heating rate Q E O is also available The magnetic field and the magnetic flux density can be plotted as vector fields using arrows Example The example demonstrates the skin effect when AC current is carried by a wire with circular cross section The conductivity of copper is 57 10 and the permeability is 1 i e y 4210 7 At the line frequency 50 Hz the w e term is negligible Due to the induction the current density in the interior of the conductor is smaller than at the outer surface where it is set to s 1 a Dirichlet condition for the electric field E 1 o For this case an analytical solution is available 2 52 Application Modes Jy kr JA J KR from the centerline and J y x is the first Bessel function of zeroth order wherek yJjouo R is the radius of the wire r is the distance Using the Graphical User Interface Start the pdet ool GUI and select the application mode AC Power Electromagnetics Draw a circle with radius 0 1 to represent a cross section of the conductor and proceed to the boundary mode to define the boundary condition Use the Select All option to select all boundaries and enter 1 57E6 intother edit field in the Boundary Condition dialog box to define the Dirichlet boundary condition E J 0 Open th
160. ively b describes the boundary conditions of the PDE problem b can be either a Boundary Condition matrix or the name of a Boundary M file The formats of the Boundary Condition matrix and Boundary M file are described in the entries on assemb and pdebound respectively adaptmesh The adapted triangular mesh of the PDE problem is given by the mesh data p e andt Details on the mesh data representation can be found in the entry on initmesh The coefficientsc a andf of the PDE problem can be given in a wide variety of ways In the context of adapt mesh the coefficients can depend on u if the nonlinear solver is enabled using the property non i n The coefficients cannot depend ont thetime A complete listing of all options can be found in the entry onassempde The table below lists the property name property value pairs their default values and descriptions of the properties Property Property Default Description Maxt positive integer inf Maximum number of new triangles Ngen positive integer 10 Maximum number of triangle generations Mesh pl el tl initmesh Initial mesh Tri pick MATLAB function pdeadworst Triangle selection method Par numeric 0 5 Function parameter Rmet hod longest regul ar longest Triangle refinement method Nonlin on of f of f Use nonlinear solver Tol n numeric le 4 Nonlinear tolerance Init u0 0 Nonlinear initial value Jac fixed lumped full fixed Nonlinear J acobian calculation nor
161. ked for correctness It is assumed that no circles or ellipses are identical or degenerated and that no lines have zero length Polygons must not be selfintersecting Use the function cs gchk tocheck the CSG model NaN is returned if the set formula sf cannot be evaluated pdetool pdegeom pdebound wgeom whound pdecirc pderect pdepoly pdeellip csgchk csgdel dst idst Purpose Synopsis Description See Also Discrete sine transform y dst x x idst y y dst x computes the discrete sine transform of the columns of x For best performance speed the number of rows in x should be 2 1 for some integer m y dst x n pads or truncates the vector x to length n before transforming Ifx isa matrix thedst operation is applied to each column x i dst y calculates the inverse discrete sine transform of the columns of y For best performance speed the number of rows in y should be 2 1 for some integer m x i dst y n pads or truncates the vector y to length n before transforming Ify is a matrix thei dst operation is applied to each column poisolv poiasma poiindex 5 35 hyperbolic Purpose Synopsis Description 5 36 Solve hyperbolic PDE problem ul hyperbolic u0 ut0 tlist b p e t c a f d ul hyperbolic u0 ut0 tlist b p e t c a f d rtol ul hyperbolic u0 ut0 tlist b p e t c a f d rtol atol ul hyperbolic u0 ut0 tlist Kk F B ud M ul hyperbolic u0
162. kg C and the thermal conductivity is 40 W m C The heat source is 20000 W m The temperature at the right end is 100 C The surrounding temperature at the outer boundary is 100 C and the heat transfer coefficient is 50 W m C The heat flux at the left end is 5000 W m But this is a cylindrical problem so you need to transform the equation using the cylindrical coordinates r z and 6 Due to symmetry the solution is independent of 6 so the transformed equation is du Of du Of du _ rp Cs Sert He fr The boundary conditions are e a kVu 5000 at the left end of the rod Neumann condition Since the generalized Neumann condition in the PDE Toolbox isa cVu qu g and c depends on r in this problem c kr this boundary condition is expressed as a cVu 5000r e u 100 at the right end of the rod Dirichlet condition e a kVu 50 100 u at the outer boundary generalized Neumann condition In the PDE toolbox this must be expressed as a cVu 50r u 50r 100 e Thecylinder axis r 0 is not a boundary in the original problem but in our 2 D treatment it has become one We must give the artificial boundary condition a cVu 0 here The initial value is u t 0 2 21 2 Examples Using the Graphical User Interface Solve this problem using the pdet ool GUI Model the rod as a rectangle with its base along the x axis and let the x axis be the z direction and the y axis be ther directi
163. l can be replaced by the boundary condition cVu Vv auv dx E qu g vds fv dx Q dQ Q The Elliptic Equation Replace the original problem with Find u such that evo Vv auv fv dx qu g vds 0 Vv This equation is called the variational or weak form of the differential equation Obviously any solution of the differential equation is alsoa solution of the variational problem The reverse is true under some restrictions on the domain and on the coefficient functions The solution of the variational problem is also called the weak solution of the differential equation The solution u and the test functions v belong to some function space V The next step is to choose an Np dimensional subspace Vy lt V Project the weak form of thedifferential equation ontoa finite dimensional function spacesimply means requesting u and vtolieinV rather than V The solution of the finite dimensional problem turns out to be the element of V that lies closest to the weak solution when measured in the energy norm see bel ow Convergence is guaranteed if the space Vy tends toV as Np gt co Since the differential operator is linear we demand that the variational equation is satisfied for N y test functions 6 e Vy that form a basis i e cVu Vo auo fojdx qu g ods 0 i 1 N Q dQ Expand u in the same basis of V y N u x Y Ujo j l and obtain the system of equations N y eve Vo
164. l to the equipotential lines of the magnetostatic potential wa PDE Toolbox MAGNET M File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help A Y Magnetostatics 0 4185 Y 0 828 Set formula C1402 C3 C4 C54R 14R2 501 Contour A Vector field B T T T ol 1 O02 i 0 2 l Info Enter axes limits Exit l AAA _ _ _ _ _ _0 _ _ _ _ _ __ _ _ __ _______ Equipotential lines and magnetic flux in two pole motor Application Modes AC Power Electromagnetics AC power electromagnetics problems are found when studying motors transformers and conductors carrying alternating currents Let s start by considering a homogeneous dielectric with coefficient of dielectricity e and magnetic permeability u with no charges at any point The fields must satisfy a special set of the general Maxwell s equations see 2 _ 2H Mor VxH Ey VxE In the absence of current we can eliminate H from the first set and E from the second set and see that both fields satisfy wave equations with wave speed Jen AE yE 0 AH ep2 H ar ar We move on to studying a charge free homogeneous dielectric with coefficient of dielectrics e magnetic permeability u and conductivity o The current density then is J oE and the waves are damped by the Ohmic resistanc
165. lbox M enus Axes Limits 1 5 1 5 Auto In the Axes Limits dialog box the range of the x axis and the y axis can be adjusted The axis range should be entered as a 1 by 2 MATLAB vector such as 10 10 If you check the Auto checkbox automatic scaling of the axis is used Pressing the Apply button applies the entered axis ranges pressing the Close button ends the Axes Limits dialog Application af From the Application submenu you can select from 10 available application modes The application modes can also be selected using the pop up menu in the upper right corner of the GUI 3 11 3 The Graphical User Interface 3 12 The available application modes are Generic Scalar the default mode Generic System Structural Mechanics Plane Stress Structural Mechanics Plane Strain Electrostatics Magnetostatics AC Power Electromagnetics Conductive Media DC Heat Transfer Diffusion See Application Modes on page 2 35 for more details PDE Toolbox M enus Draw Menu Draw Boundary PDE Mesh Solve Plot Window Draw Mode Rectangle square Rectangle square centered Ellipse circle Ellipse circle centered Polygon Rotate Export Geometry Description Set Formula Labels Draw Mode Enter draw mode Rectangle square Draw a rectangle square starting at a corner Using the left mouse button click and drag to create a rectangle Using the right mouse button or Control
166. les of Elliptic Problems Also initialize a vector error for the maximum norm errors of the successive solutions and set the initial error err to1 Theloop then runs until the error of the solution is smaller than 10 1 Refinethe mesh The current triangular mesh defined by the geometry circl eg the point matrix p the edge matrix e and thetriangle matrixt is refined and the mesh is returned using the same matrix variables 2 Assemble and solve the linear system Note that the coefficients of the elliptic PDE are constants c f 1 a 0 for this simple case ci rcleb1 contains a description of the boundary conditions andp e andt definethe triangular mesh 3 Find the error of the numerical solution produced by the PDE Toolbox The vector exact contains the exact solution at the nodes Note that what you actually find is the max norm error of the solution at the nodes 4 Plot the mesh the solution and the error Notice that the plot function pdesurf asthird argument can take any vector of values on the mesh given byp andt not just the solution In this case you are also plotting the error function x10 7 SN N NS AW A INN AN YY X N LO LA A A f M iy NS NNSSS SH NS IS N Y Y AN Y N SS IRN RY gt IRIS ASSA NORRIS 0 SR VAVA SS A 0 5 N NY A NY y INN S Ny A A VAVA Y a A A AVA N N N NY S ZA As Ly A Al Ai A ZA
167. ll subdomain borders by selecting Remove All Subdomain Borders fromthe Boundary menu Thetwo boundaries at the inset in the lower left are clamped i e Dirichlet conditions with zero displacements The rounded cut is subject toa Neumann condition with q 0 andg1 0 5 nx g2 0 5 ny The remaining boundaries are free no normal stress that is a Neumann condition with q 0 andg 0 The next step is to open the PDE Specification dialog box and enter the PDE parameters ka PDE Specification Equation Structural mechanics plane stress Type of PDE Coefficient Value Description Elliptic E 196E3 Young s modulus P nu 0 31 Poisson ratio Kx 0 Volume force x direction wv Eigenmodes Ky 0 Volume force y direction rho 1 0 Density OK Cancel al TheE and v nu parameters are Young s modulus and Poisson s ratio respectively There are no volume forces so Kx andKy are zero p rho is not used in this mode The material is homogeneous sothe same E and v apply to the whole 2 D domain Initialize the mesh by pressing the A button If you want you can refine the mesh by pressing the Refine button The problem can now be solved by pressing the button A number of different strain and stress properties can be visualized such as the displacements u and v the x and y direction strains and stresses the shear stress the von Mises effective stress and the principal stresses and strains All t
168. ller than 100 Plot the first eigenmode and compare it to MATLAB s membrane function gt pdesurf p t v 1 figure membrane 1 20 9 9 2 28 Examples of Eigenvalue Problems membrane can produce the first 12 eigenfunctions for the L shaped membrane Compare also the 12th eigenmodes figure gt pdesurf p t v 12 figure gt membrane 12 20 9 9 Looking at the following eigenmodes you can see how the number of oscillations increases The eigenfunctions are symmetric or antisymmetric around the diagonal from 0 0 to 1 1 which divides the L shaped membrane into two mirror images In a practical computation you could take advantage of such symmetries in the PDE problem and solve over a region half the size The eigenvalues of the full L shaped membrane are the union of those of the half with Dirichlet boundary condition along the diagonal eigenvalues 2 4 7 11 13 16 and 17 and those with Neumann boundary condition eigenvalues 1 3 5 6 10 12 14 and 15 The eigenvalues Ag and Ag make up a double eigenvalue for the PDE at around 49 64 Also the eigenvalues 4 2 and A make up another double eigenvalue at around 99 87 You may have gotten two different but close values The default triangulation made by i ni t mesh is not symmetric around the diagonal but a symmetric grid gives a matrix with a true double eigenvalue Each of the eigenfunctions ug and uy consists of three copies of eigenfunction
169. lly started and the polygon is drawn in an empty geometry model The optional argument abel assigns a name to the polygon otherwise a default name is chosen The state of the Geometry Description matrix inside pdet oo is updated to include the polygon You can export the Geometry Description matrix from pdetool by using the Export Geometry Description option from the Draw menu The format of the Geometry Description matrix is described in the entry on decsg The command pdepoly 1 0 0 1 1 1 0 0 1 1 1 1 creates the L shaped membrane geometry as one polygon pdeci rc pderect pdet ool 5 71 pdeprtni Purpose Synopsis Description Cautionary See Also 5 72 I nterpolate from triangle midpoint data to node data un pdeprtni p t ut un pdeprtni p t ut gives linearly interpolated values at node points from the values at triangle midpoints The geometry of the PDE problem is given by the mesh datap andt Details on the mesh data representation can be found in the entry oni nit mesh Let N be the dimension of the PDE system np the number of node points and nt the number of triangles The components of triangle data inut are stored as N rows of length n The components of the node data are stored in un as N columns of length np pdeprtni andpdei ntrp arenot inverse functions Theinterpolation introduces some averaging assempde initmesh pdeintrp pderect Purpose Synopsis Descr
170. lux vector field plots Example In the following example a heat transfer problem with differing material parameters is solved The problem s 2 D domain consists of a square with an embedded diamond a square with 45 degrees rotation The square region consists of a material with coefficient of heat conduction of 10 and a density of 2 The diamond shaped region contains a uniform heat source of 4 and it has a coefficient of heat conduction of 2 and a density of 1 Both regions have a heat capacity of 0 1 Application Modes Using the Graphical User Interface Start the pdet ool GUI and select the application mode Heat Transfer n draw mode set the x and y axis limits to 0 5 3 5 and select the Axis Equal option from the Options menu The square region has corners in 0 0 3 0 3 3 and 0 3 The diamond shaped region has corners in 1 5 0 5 2 5 1 5 1 5 2 5 and 0 5 1 5 The temperature is kept at 0 on all the outer boundaries so you do not have to change the default boundary conditions M ove on to define thePDE parameters make sure that the Heat Transfer application mode is selected in the PDE mode by double clicking on each of the two regions and enter the PDE parameters You want to solve the parabolic heat equation so make sure that the Parabolic option is selected In the square region enter a density of 2 a heat capacity of 0 1 and a coefficient of heat conduction of 10 There is no heat source so set
171. lve Plot Window Help OF ti X 0 4186 Y 0 828 Setformula C1 C2 C3 C44C5 A1 R2 501 pl L L 1 L 1 L n 02 01 0 0 02 03 04 05 06 07 08 09 1 11 12 Info Click and drag at perimeter to create ellipse Exit Before moving tothe PDE mode select the boundaries along the x axis and set the boundary condition to a Neumann condition with g 0 and q 0 In the PDE mode turn on the labels by selecting the Show Subdomain Labels option from the PDE menu Double click on each subdomain to definethePDE parameters e Inthecoil both u andj are1 sothe default values do not need to be changed e In the stator and the rotor u is nonlinear and defined by the equation above Enter u as 5000 1 0 05 ux 2 uy 2 200 ux 24uy 2 is equal to VA is O no current e Intheair gap uis 1 and is 0 2 49 2 Examples I nitializethe mesh and continue by opening the Solve Parameters dialog box by selecting Parameters from the Solve menu Since this is a nonlinear problem thenonlinear solver must beinvoked by checking the Use nonlinear solver If you want you can adjust the tolerance parameter The adaptive solver can be used together with the nonlinear solver Solve the PDE and plot the magnetic flux density B using arrows and the equipotential lines of the magnetostatic potential A using a contour plot The plot clearly shows as expected that the magnetic flux is paralle
172. m norm NOTE The adaptive mode and the nonlinear solver can be used together e Parabolic PDEs The solve parameters for the parabolic PDEs are Time A MATLAB vector of times at which a solution to the parabolic PDE should be generated Therelevant timespan is dependent on the dynamics of the problem Examples 0 10 and ogspace 2 0 20 u t0 The initial value u t for the parabolic PDE problem The initial value can be a constant or a column vector of values on the nodes of the current mesh Relative tolerance Relative tolerance parameter for the ODE solver that is used for solving the time dependent part of the parabolic PDE problem Absolute tolerance Absolute tolerance parameter for the ODE solver that is used for solving the time dependent part of the parabolic PDE problem 3 27 3 The Graphical User Interface 3 28 Solve Parameters gt u tO 01 bsolute tolerance lative tolerance 001 OK Cancel l Solve Parameters dialog box for hyperbolic PDEs e Hyperbolic PDEs The solve parameters for the hyperbolic PDEs are Time A MATLAB vector of times at which a solution tothe hyperbolic PDE should be generated The relevant time span is dependent on the dynamics of the problem Examples 0 10 andlogspace 2 0 20 u t0 The initial value u t for the hyperbolic PDE problem The initial value can be a constant or a column vec
173. m numerid inf energy inf Nonlinear residual norm Par is passed totheTri pick function TheTri pick function is described below Normally it is used as tolerance of how well the solution fits the equation 5 9 adaptmesh 5 10 No morethan Ngen successive refinements are attempted Refinement is also stopped when the number of triangles in the mesh exceeds Maxt p1 el andt1 aretheinput mesh data This triangular mesh is used as starting mesh for the adaptive algorithm Details on the mesh data representation can be found in the entry oni nit mesh If noinitial mesh is provided the result of a Call toi nit mesh with no options is used as initial mesh The triangle selection method Tri pi ck is a user definable triangle selection method Given the error estimate computed by the function pdej mps the triangle selection method selects the triangles to be refined in the next triangle generation The function is called usingtheargumentsp t cc aa ff u errf andpar p andt represent the current generation of triangles cc aa andf f are the current coefficients for the PDE problem expanded to triangle midpoints u is the current solution errf is the computed error estimate and par the function parameter given toadapt mesh as optional argument The matricescc aa ff anderrf all have N columns where N is the current number of triangles The number of rowsincc aa andff are exactly the same as the input arguments c a andf errf ha
174. m the pop up menu in the upper right part of the GUI or by selecting an application from the Application submenu in the Options menu Note that changing the application resets all PDE coefficients and boundary conditions to the default values for that specific application mode 2 35 2 Examples When using an application mode the generic PDE coefficients are replaced by application specific parameters such as Y oung s modulus for problems in structural mechanics The application specific parameters are entered by selecting Parameters fromthe PDE menu or by pressing the PDE button You can also access the PDE parameters by double clicking on a subdomain if you arein the PDE mode That way it is possible to define PDE parameters for problems with regions of different material properties The Boundary condition dialog box is also altered so that the Description column reflects the physical meaning of the different boundary condition coefficients Finally the Plot Selection dialog box allows you to visualize the relevant physical variables for the selected application NOTE In the User entry options in the Plot Selection dialog box the solution and its derivatives are always referred toasu ux anduy v vx and vy for the system cases even if the application mode is nongeneric and the solution of the application specific PDE normally is named e g V or T In the remaining part of this section we explain each of the application m
175. mains remove borders between subdomains and export the decomposed geometry and the boundary conditions to the workspace PDE menu The PDE menu provides a dialog box for specifying the PDE and there are menu options for labeling subdomains and exporting PDE coefficients to the workspace 1 24 Using the Graphical User Interface Mesh menu From the Mesh menu you create and modify the triangular mesh Y ou can initialize refine and jiggle the mesh undo previous mesh changes label nodes and triangles display the mesh quality and export the mesh to the workspace Solve menu From the Solve menu you solve the PDE You can also open a dialog box where you can adjust the solve parameters and you can export the solution to the workspace Plot menu From the Plot menu you can plot a solution property A dialog box lets you select which property to plot which plot style to use and several other plot parameters If you have recorded a movie animation of the solution you can export it to the workspace Window menu The Window menu lets you select any currently open MATLAB figure window The selected window is brought to the front Help menu The Help menu provides a brief help window The Toolbar The toolbar underneath the main menu at the top of the GUI contains icon buttons that provide quick and easy access to some of the most important functions O B Se gt a0 ref A a lO The five left most buttons are Draw mode
176. n angles By default the rotation center is the center of mass of the selected objects If the Use center of mass option is unchecked you can enter a rotation center Xe Yc as a 1 by 2 MATLAB vector such as 0 4 0 3 PDE Toolbox M enus Boundary Menu Boundary PDE Mesh Solve Plot Window Help Boundary Mode lt Ctrt gt b Specify Boundary Conditions Show Edge Labels Show Subdomain Labels 3 AH Gul ders Export D ecomposed Geometry Boundary Cond s Boundary Mode Specify Boundary Conditions Show Edge Labels Show Subdomains Labels Remove Subdomain Border Remove All Subdomain Borders Export Decomposed Geometry Boundary Cond s Enter the boundary mode Specify boundary conditions for the selected boundaries If no boundaries are selected the entered boundary condition applies to all boundaries Toggle the labeling of the edges outer boundaries and subdomain borders on off The edges are labeled using the column number in the Decomposed Geometry matrix Toggle the labeling of the subdomains on off The subdomains are labeled using the subdomain numbering in the Decomposed Geometry matrix Remove selected subdomain borders Remove all subdomain borders Export the Decomposed Geometry matrix g and the Boundary Condition matrixb tothe main workspace 3 15 3 The Graphical User Interface 3 16 Specify Boundary Conditions al Specify b
177. n depend ont the time A complete listing of all options can be found in the entry onassempde hyperbolic Examples See Also atol andrtol areabsoluteand relative tolerances that are passed to the ODE sol ver ul hyperbolic u0 ut0 tlist K F B ud M produces the solution tothe ODE problem 2 d u B MB K u F u Bu uy dt with initial values for u beingu0 andut0 Solve the wave equation 2 du ar Au on a square geometry 1 lt x y lt 1 s quar eg with Dirichlet boundary conditions u 0 for x 1 and Neumann boundary conditions du F for y 1 s quar eb3 Choose u 0 arctan cos 1x and Lo 3sin zx exp cos ry Compute the solution at times 0 1 6 1 3 29 6 5 t initmesh squareg p ane p ae O atan cos pi 2 x utO 3 sin pi x explcos pi y tlist linspace 0 5 31 uu hyperbolic u0 ut0 tlist squareb3 p e t 1 0 0 1 p e p 1 p 2 X y u The file pdedemo6 contains a complete example with animation assempde parabolic 5 37 initmesh Purpose Synopsis Description 5 38 Create initial triangular mesh p e t i nit mesh g p e t initmesh g PropertyName PropertyValue p e t initmesh g returns a triangular mesh using the geometry specification function g It uses a Delaunay triangulation algorithm The mesh size is determined from the shape of the geometry y describes the geometry of the
178. n matrix gd are valid given optional real numbers x i mandy i mas current length of the x and y axis and using a special format for polygons F or a polygon the last vertex coordinate can be equal to thefirst one toindicatea closed polygon If x i mandy i m are specified and if the first and the last vertices are not equal the polygon is considered as closed if these vertices are within a certain closing distance These optional input arguments are meant to be used only when callingcsgchk frompdetool gstat csgchk gd is identical to the above call except for using the same format of gd that is used by decsg This call is recommended when using csgchk aS a command line function gstat iS a row vector of integers that indicates the validity status of the corresponding solid objects i e columns in gd For a circle solid gst at 0 indicates that the circle has a positive radius 1 indicates a nonpositive radius and 2 indicates that the circle is not unique For a polygon gst at O indicates that the polygon is closed and does not intersect itself i e it has a well defined unique interior region 1 indicates an open and non self intersecting polygon 2 indicates a closed and self intersecting polygon and 3 indicates an open and self intersecting polygon For a rectangle solid gst at is identical to that of a polygon This is so because a rectangle is considered as a polygon by csgchk For an ellipse solid gst at 0
179. n the set formula the name refers to the set of points inside the object The resulting geometrical model is the set of points for which the set pdetool formula evaluates to true For a description of the syntax of the set formula seedecsg By default the resulting geometrical model is the union of all objects A snap to grid function is available This means that objects align tothe grid The grid can be turned on and off and the scaling and the grid spacing can be changed In boundary mode you can specify the boundary conditions You can have different types of boundary conditions on different boundaries In this mode the original shapes of the solid building objects constitute borders between subdomains of the model Such borders can be eliminated in this mode The outer boundaries are color coded to indicate the type of boundary conditions A red outer boundary corresponds to Dirichlet boundary conditions blue to generalized Neumann boundary conditions and green to mixed boundary conditions You can return tothe boundary condition display by clicking on the dQ button or by selecting Boundary Mode from the Boundary menu In PDE mode you can specify the type of PDE problem and the coefficients c a f and d You may specify the coefficients for each subdomain independently This makes it easy to specify e g various material properties in one PDE model The PDE to be solved can be specified by clicking on the PDE button
180. nberg Stenger 8 is used in which the longest side of atriangleis always split whenever any of the sides have hanging nodes This guarantees that no angleis ever smaller than half the smallest angle of the original triangulation Two selection criteria can be used One pdeadworst refines all elements with value of the error indicator larger than half the worst of any element The other pdeadgsc refines all elements with an indicator value exceeding a user defined dimensionless tolerance The comparison with the tolerance is properly scaled with respect to domain and solution size etc The Termination Criteria For smooth solutions error equidistribution can be achieved by thepdeadgsc selection if the maximum number of elements is large enough Thepdeadworst adaption only terminates when the maximum number of elements has been exceeded This mode is natural when the solution exhibits singularities The error indicator of the elements next to the singularity may never vanish regardless of element size and equidistribution is too much to hope for 4 28 Fast Solution of Poisson s Equation Fast Solution of Poisson s Equation While the general strategy of the PDE Toolbox is to use MATLAB s built in solvers for sparse systems there are situations where faster solution algorithms are available Onesuch exampleis when solving Poisson s equation Au f in Q with Dirichlet boundary conditions where Q is a rectangle For the fast
181. nd P P Gis oP eal i 1 2 where Pp is the mid point of P P 5 The Elliptic Equation For each triangle the vertices Pm of the local triangle correspond to the indices im Of the mesh points The contributions of the individual triangle are added to the matrices such that e g K tK ln ict Rene m n 1 2 3 This is done by the function ass empde The gradients and the areas of the triangles are computed by the function pdetrg The Dirichlet boundary conditions are treated in a slightly different manner They are eliminated from the linear system by a procedure that yields a symmetric reduced system The function assempde can return matrices K F B and ud such that the solution is u Bv ud where Kv F u is an Np vector and if therank of the Dirichlet conditions is rD then v has Np rD components 4 9 4 The Finite Element M ethod The Elliptic System Thetoolbox can also handle systems of N partial differential equations over the domain Q We have the elliptic system V c Vu au f the parabolic system qn V c Vu au f the hyperbolic system du d V c Vu au f Of and the eigenvalue system V c Vu au Adu where e is an N by N by 2 by 2 tensor By the notation V c Vu we mean the N by 1 matrix with i 1 component Fan as E E Nox Lj L lox Ax 129y dy 42 lax ay bh 2ay Jd J The symbols a and d denote N by N matrices and f and u denote column vectors of
182. ng the PDE Thenew solution however can be plotted using this dialog box PDE Toolbox M enus F or the parabolic and hyperbolic PDEs the bottom right part of the dialog box contains an additional option Time for plot A pop up menu allows you to select which of the solutions to plot by selecting the corresponding time By default the last solution is plotted 3 35 3 The Graphical User Interface 3 36 Also the Animation plot type is enabled In its property field you find an Options push button If you press on it an additional dialog box appears It contains parameters that control the animation Animation rate fps For the animation this parameter controls the speed of the movie in frames per second fps Number of repeats The number of times the movie is played Replay movie If you check this option the current movie is replayed without re recording the movie frames If there is no current movie this option is disabled vj Plot Selection Plot type Property User entry Plot style W Color r patacas E N interpolated shad Contour Arrows grad u za proportional 4 I Deformed mesh grad u Height 3 D plot u continuous 4 J A Dp Plot in x y grid Contour plot levels 20 Ill Plot solution automatically 1 Show mesh Colormap cool Eigenvalue 9966 Piot Done Canel al For eigenvalue problems the bottom ri
183. nge The default entry is 0 100 which is just what you want 2 27 2 Examples Finally solve the L shaped membrane problem by pressing the button The solution displayed is the first eigenfunction The value of the first smallest eigenvalue is also displayed You find the number of eigenvalues on the information line at the bottom of the GUI You can open the Plot Selection dialog box and choose which eigenfunction to plot by selecting from a pop up menu of the corresponding eigenvalues Using Command Line Functions The geometry of the L shaped membrane is described in the file shapeg mand the boundary conditions in the file shapeb m First initialize the mesh and refine it twice using the command line functions at the MATLAB prompt p e t initmesh Ishapeg p e t refinemesh shapeg p e t p e t refinemesh shapeg p e t Recall the general eigenvalue PDE problem description V cVu au Adu This means that in this case you have c 1 a 0 and d 1 The syntax of pdeei g the eigenvalue solver in the PDE Toolbox is v pdeeig b p e t c a d r The input argument r is a two element vector indicating the interval on the real axis wherepdeeig searches for eigenvalues Here you are looking for eigenvalues lt 100 so the interval you useis 0 100 Now you can call pdeei g and see how many eigenvalues you find v 1 pdeeig lshapeb p e t 1 0 1 0 100 There are 19 eigenvalues sma
184. nl n2 K poi as ma n K poiasma n1 n2 h1 h2 assembles the contributions to the stiffness matrix from boundary points n1 and n2 are the numbers of points in the first and second directions andh1 and h2 are the mesh spacings K is a sparse n1 n2 by n1 n2 matrix The point numbering is the canonical numbering for a rectangular mesh K poiasma nl n2 useshl h2 K poiasma n USeSnl n2 n poiindex poisolv 5 85 poicalc Purpose Synopsis Description See Also 5 86 Fast solver for Poisson s equation on a rectangular grid u poicalc f hl h2 nl n2 u poicalc f hl h2 u poicalc f u poicalc f h1 h2 n1 n2 calculates the solution of Poisson s equation for the interior points of an evenly spaced rectangular grid The columns of u contain the solutions corresponding to the columns of the right hand sidef h1 andh2 arethe spacings in thefirst and second direction andn1 andn2 arethe number of points The number of rows inf must ben1 n2 Ifn1 andn2 are not given the square root of the number of rows off is assumed Ifh1 andh2 are not given they are assumed to be equal The ordering of the rowsinu andf is the canonical ordering of interior points as returned by poi index The solution is obtained by sine transforms in the first direction and tridiagonal matrix solution in the second direction n1 should be 1 less than a power of 2 for best performance poiindex poiasma dst idst poisolv polindex Purpose
185. ns Using Command Line Functions Although thepdet 00 GUI provides a convenient working environment there are situations where the flexibility of using the command line functions is needed These include e Geometrical shapes other than straight lines circular arcs and elliptical arcs Nonstandard boundary conditions Complicated PDE or boundary condition coefficients e More than two dependent variables in the system case e Nonlocal solution constraints e Special solution data processing and presentation itemize The GUI can still be a valuable aid in some of the situations presented above if part of the modeling is done using the GUI and then made available for command line use through the extensive data export facilities of the GUI Data Structures and Utility Functions The process of defining your problem and solving it is reflected in the design of the GUI A number of data structures define different aspects of the problem and the various processing stages produce new data structures out of old ones See the figure on the next page The rectangles are functions and ellipses are data represented by matrices or M files Arrows indicate data necessary for the functions As thereis a definite direction in this diagram you can cut intoit by presenting the needed data sets and then continue downward In the sections below we give pointers to descriptions of the precise formats of the various data structures and M files 1
186. ns and Neumann boundary conditions are also called natural boundary conditions See Chapter 4 TheF inite Element Method for the general system case 1 4 Intro duction In Which Areas Can the Toolbox Be Used The PDEs implemented in the toolbox are used as a mathematical model for a wide variety of phenomena in all branches of engineering and science The following is by no means a complete list of examples The elliptic and parabolic equations are used for modeling e steady and unsteady heat transfer in solids e flows in porous media and diffusion problems e electrostatics of dielectric and conductive media e potential flow The hyperbolic equation is used for e transient and harmonic wave propagation in acoustics and electromagnetics e transverse motions of membranes The eigenvalue problems are used for e g e determining natural vibration states in membranes and structural mechanics problems Last but not least the toolbox can be used for educational purposes as a complement to understanding the theory of the Finite Element Method How Do I Define a PDE Problem Thesimplest way to definea PDE problem is using the graphical user interface GUI implemented in pdet ool There are three modes that correspond to different stages of defining a PDE problem e Draw mode you create the geometry using the constructive solid geometry CSG model paradigm A set of solid objects rectangle circle ellipse and
187. ns can only be activated one at the time and they all work the same way single clicking on a button allows you to draw one solid object of the selected type Double clicking on a button makes it stick and you can then continue to draw solid objects of the selected type until you single click on the button to release it If you are not in the draw mode when you click on one of the draw mode buttons the GUI enters the draw mode automatically The Toolbar The second group of buttons includes the following buttons from left to right e JQ Enters the boundary mode e PDE Opens the PDE Specification dialog box e A Initializes the triangular mesh e The Refine button the four triangle icon Refines the triangular mesh e Solves the PDE e The 3 D solution plot icon Opens the Solution Plot Selection dialog box The buttons in the second group are of the flash type single clicking on them initiates the associated function The last right most button with the magnifier icon toggles the zoom function on off 3 39 3 The Graphical User Interface 3 40 The Finite Element Method The Elliptic Equation 0 2 eee ee eee 4 3 The Elliptic System 0 0000 4 10 The Parabolic Equation 0008 4 13 The Hyperbolic Equation 0 00 4 16 The Eigenvalue Equation 048 4 17 Nonlinear Equations
188. nsionality reduction When the problem is such that variation with z is negligible all z derivatives drop out and the 2 D equation has exactly the same units and coefficients as in 3 D Slab geometries are treated by integration through the thickness The result is a 2 D equation for the z averaged solution with the thickness say D x y multiplied onto all the PDE coefficients c a d and f etc For instance if you want to compute the stresses in a sheet welded together from plates of different thickness multiply Young s modulus E volume forces and specified surface tractions by D x y Similar definitions of the equation coefficients are called for in other slab geometry examples and application modes 1 8 G etting Started Getting Started To get you started let s use the graphical user interface GUI pdet oo which isa part of the PDE Toolbox tosolve a PDE step by step The problem that we would like to solve is Poisson s equation Au f The 2 D geometry on which we would like to solve the PDE is quite complex The boundary conditions are of Dirichlet and Neumann types First invoke MATLAB To start the GUI type the command pdet 001 at the MATLAB prompt It can take a minute or two for the GUI to start The GUI looks similar to the figure below with exception of the grid Turn on the grid by selecting Grid from the Options menu Also enable the snap to grid feature by selecting Snap from the Options menu
189. nt pdesurf 5 65 pdenonlin Purpose Synopsis Description 5 66 Solve nonlinear PDE problem u res pdenonlin b p e t c a f u res pdenonlin b p e t c a f PropertyName PropertyValue u res pdenonlin b p e t c a f solves the nonlinear PDE scalar PDE problem V cVu au f on Q or the nonlinear system PDE problem V c 8 Vu au f on 9 where the coefficients c a and f may depend on u The algorithm solves the equation by using damped Newton iteration with the Armijo Goldstein line search strategy The solution u is represented as the solution vector u Details on the representation of the solution vector can befoundintheentry onassempde res is the norm of the Newton step residuals The triangular mesh of the PDE problem is given by the mesh data p e andt Details on the mesh data representation can be found in the entry oni ni t mesh b describes the boundary conditions of the PDE problem b can be either a Boundary Condition matrix or the name of a Boundary M file The formats of the Boundary Condition matrix and Boundary M file are described in the entries on assemb and pdebound respectively respectively Note that for the general call to pdebound the boundary conditions can also depend on u A fixed point iteration strategy is employed to solve for the nonlinear boundary conditions The coefficients c a f of the PDE problem can be given in a wide variety of ways In the contex
190. ntains du Z dy There is one column for each triangleint in both ux anduy The geometry of the PDE problem is given by the mesh data p andt Details on the mesh data representation can be found in the entry oni nit mesh The format for the solution vector u is described inassempde The optional argument sd restricts the computation to the subdomains in the lists dl assempde 5 61 pdeintrp Purpose Synopsis Description Cautionary See Also 5 62 Interpolate from node data to triangle midpoint data ut pdeintrp p t un ut pdeintrp p t un gives linearly interpolated values at triangle midpoints from the values at node points The geometry of the PDE problem is given by the mesh datap andt Details on the mesh data representation can be found in the entry oni nit mesh Let N be the dimension of the PDE system np the number of node points and n the number of triangles The components of the node data are stored in un either as N columns of length n or as an ordinary solution vector The first np values of un describe the first component the following np values of un describe the second component and so on The components of triangle data are stored inut as N rows of length nj pdeprtni andpdei ntrp arenot inverse functions Theinterpolation introduces some averaging assempde initmesh pdeprtni pdejmps Purpose Synopsis Description See Also Error estimates for adaption errf pdejmps p
191. nting the circular conductors Place them symmetrically with centers in 0 6 0 and 0 6 0 If the rectangle s label is R1 and the circles labels are C1 and C2 the 2 D domain of the problem is expressed by the set formula R1 C1 C2 Enter the set formula and press the dQ button to decompose the geometry and enter the boundary mode Select all the outer boundaries and enter the Neumann boundary condition into the Boundary Condition dialog box For the left circular conductor boundaries enter the Dirichlet boundary condition V 1 and for the right circular conductor enter the Dirichlet condition V 1 Next open the PDE Specification dialog box and enter 0 into the edit box for the current source q The default value for the conductivity o is 1 and needs not to be changed Initialize the mesh and refine it twice End by jiggling the mesh once to improve the triangle quality Solve the PDE by pressing the button The resulting potential is zero along the y axis which is a vertical line of anti symmetry for this problem 2 56 Application Modes Visualize the current density J by plotting the absolute value using a contour plot and the vector field using arrows The current flows as expected from the conductor with a positive potential tothe conductor with a negative potential Y PDE Toolbox untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help N Conductive Media DC X 0 8
192. o the re entrant corners where the gradients are steep For example use the triangle selection method picking the worst triangles and set the maximum number of triangles to 500 Add one uniform mesh refinement by pressing the Refine button once Finally turn adaptive mode off and press the button once more 2 44 Application Modes Tolook at the equipotential lines select a contour plot from the Plot Selection dialog box To display equipotential lines at every 100th volt enter 0 100 1000 into the Contour plot levels edit box xf PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help 0 09222 Y 0 262 0 25 A L L 1 1 1 1 1 1 1 1 1 L i L 0 4 0 35 0 3 0 25 0 2 0 15 0 1 0 05 0 0 05 01 015 02 025 0 3 0 35 Into Select a different plot or change mode to alter PDE mesh or boundary conditions Equipotential lines in air filled frame 2 Examples Magnetostatics Magnets electric motors and transformers are areas where problems involving magnetostatics can be found The statics implies that the time rate of change is slow so we start with Maxwell s equations for steady cases VxH J v B 0 and the relationship B uH where B is the magnetic flux density H is the magnetic field intensity J is the current density and u is the material s magn
193. odes in more detail and give examples of how to solve application specific problems using the PDE Toolbox Structural Mechanics Plane Stress In structural mechanics the equations relating stress and strain arise fromthe balance of forces in the material medium Planestress is a condition that prevails in a flat platein the xy plane loaded only in its own plane and without z direction restraint The stress strain relation can then be written assuming isotropic and isothermal conditions Ox E ly 0 Ex Oe rw vil 0 E vV Te 00 1 v 2 Tiy where o and oy are the stresses in the x and y directions and ty is the shear stress The material properties are expressed as a combination of E the dastic modulus or Young s modulus andv Poisson s ratio 2 36 Application Modes The deformation of the material is described by the displacements in the x and y directions u and v from which the strains are defined as Me WV _ du av x Ox y dy ley dy dx The balance of force equations are 90 OT eS ee A ox oy z dT 00 xy Ys o G where Kx and Ky are volume forces body forces Combining the relations above wearriveat the displacement equations which can be written V c Vu k where c is a rank four tensor which can be written as four 2 by 2 matrices Cir C12 Gy and Cp E ce 1 7 0 G 0 a OH 0G cano Cy G 0 0 2G Hu where G the shear modulus is defined by E 21 v 2 37
194. of the property that is represented using color or contour lines e Property contains four pop up menus containing lists of properties that are available for plotting using the corresponding plot type From the first pop up menu you control the property that is visualized using color and or contour lines The second and third pop up menus contain vector valued properties for visualization using arrows and deformed mesh respectively From the fourth pop up menu finally you control which scalar property to visualize using z height in a 3 D plot The lists of properties are dependent on the current application mode For the generic scalar mode you can select the following scalar properties u The solution itself abs grad u The absolute value of Vu evaluated at the center of each triangle abs c grad u The absolute value of c Vu evaluated at the center of each triangle user entry A MATLAB expression returning a vector of data defined on the nodes or the triangles of the current triangular mesh The solution u its derivatives ux and uy the x and y components of c Vu cux and cuy and x and y areall available in the local workspace Y ou enter the expression into the edit box to the right of the Property pop up menu in the User entry column Examples u u x y The vector property pop up menus contain the following properties in the generic scalar case grad u The negative gradient of u Vu c grad u c times the ne
195. omain and to define boundary conditions for a PDE problem It also makes it possible to specify the partial differential equation to create inspect and refine the mesh and to compute and display the solution from the GUI pdetool contains several different modes In Draw mode you construct a Constructive Solid Geometry model CSG model of the geometry Y ou can draw solid objects that can overlap There are four types of solid objects e Circle object represents the set of points inside a circle Polygon object represents the set of points inside the polygon given by a set of line segments e Rectangle object represents the set of points inside the rectangle given by a set of line segments e Ellipse object represents the set of points inside an ellipse The ellipse can be rotated The solid objects can be moved and rotated Operations apply to groups of objects by doing multiple selects A Select all option is also available You can cut and paste among the selected objects The model can be saved and restored pdetool can be started by just typing the name of the model This starts the corresponding M file that contains the MATLAB commands necessary to create the model The solid objects can be combined by typing a set formula Each object is automatically assigned a unique name which is displayed in the GUI on the solid object itself The names refer to the object in the set formula More specifically i
196. on A rectangle with corners in 1 5 0 1 5 0 1 5 0 2 and 1 5 0 2 would then model a rod with length 3 and radius 0 2 Enter the boundary conditions by double clicking on the boundaries to open the Boundary Condition dialog box For the left end use Neumann conditions with 0 for q and5000 y for g For theright end use Dirichlet conditions with 1 forh and100 forr For the outer boundary use Neumann conditions with 50 y for q and50 y 100 for g For the axis use Neumann conditions with 0 for q andg Enter the coefficients into the PDE Specification dialog box is40 y a is zero d is 7800 500 y andf is20000 y Animate the solution over a span of 20000 seconds computing the solution every 1000 seconds We can see how heat flows in over the right and outer boundaries as long as u lt 100 and out when u gt 100 You can also open the PDE Specification dialog box and change the PDE type to Elliptic This shows the solution when u does not depend on time i e the steady state solution The profound effect of cooling on the outer boundary can be demonstrated by setting the heat transfer coefficient to zero 2 22 Examples of Hyperbolic Problems Examples of Hyperbolic Problems This section describes the solution of some hyperbolic PDE problems The problems are solved using the graphical user interface and the command line functions of the PDE Toolbox The Wave Equation As an example of a hyperbolic PDE let s solve th
197. onorcrn 2 12 Examples of Parabolic Problems 2 16 The Heat Equation A Heated Metal Block 2 16 Heat Distribution in Radioactive Rod ooooooooo 2 21 Examples of Hyperbolic Problems 2 23 The Wave Equation us enere tene 0 00 c eee tee 2 23 Examples of Eigenvalue Problems 2 27 Eigenvalues and Eigenfunctions for the L Shaped Membrane 20 00 c eee eee 2 27 L Shaped Membrane with Rounded Corner 2 31 Eigenvalues and Eigenmodes of aSquare 2 32 Application Modes 00 0 2 35 The Application Modes and the GUI 2 055 2 35 Structural Mechanics Plane Stress 00055 2 36 Structural Mechanics Plane Strain 2 05 2 41 Electrostatics tibiae dy Pet edt itd Sate Bite 2 43 MagnetostaticS 0 2 cee 2 46 AC Power Electromagnetics 020 e eee eee ee 2 51 Conductive Media DC 0 00 cee ee 2 55 Heat Transfer 2 22 02 NS 2 57 DIFUSIO can ot crane hes eer ee tee Oe NY 2 61 2 Examples Examples of Elliptic Problems This section describes the solution of some elliptic PDE problems The last problem a minimal surface problem is nonlinear and illustrates the use of the nonlinear solver The problems are solved using both the graphical user interface and the command line functions of the PDE Toolbox Poisson s Equation on Unit Disk As a first
198. ontexts the generalized Neumann boundary conditions is also referred to as the Robin boundary conditions In variational calculus Dirichlet conditions are also called essential boundary conditions and restrict the trial space Neumann conditions are also called natural conditions and arise as necessary conditions for a solution The variational form of the PDE Toolbox equation with Neumann conditions is given below The approximate solution to the elliptic PDE is found in three steps 1 Describe the geometry of the domain Q and the boundary conditions This can be done either interactively using pdet oo see Chapter 3 The Graphical User Interface or through M files see the entries on pdegeom and pdebound in Chapter 5 Reference 2 Builda triangular mesh on the domain Q The toolbox has mesh generating and mesh refining facilities A mesh is described by three matrices of fixed format that contain information about the mesh points the boundary segments and the triangles 4 3 4 The Finite Element M ethod 3 Discretizethe PDE and the boundary conditions to obtain a linear system Ku F The unknown vector u contains the values of the approximate solution at the mesh points the matrix K is assembled from the coefficients c a h and q and the right hand side F contains essentially averages of f around each mesh point and contributions from g Once the matrices K and F are assembled you havethe entire MATLAB environment at yo
199. opsis Description 5 90 Refine a triangular mesh pl el t1 refinemesh g p e t pl el tl refinemesh g p e t regular pl el tl refinemesh g p e t longest pl el tl refinemesh g p e t it pl el tl refinemesh g p e t it regular pl el tl refinemesh g p e t it longest pl el tl ul refinemesh g p e t u pl el tl ul refinemesh g p e t u regular pl el tl ul refinemesh g p e t u longest pl el tl ul refinemesh g p e t u it pl el tl ul refinemesh g p e t u it regular pl el tl ul refinemesh g p e t u it longest pl el tl refinemesh g p e t returns a refined version of the triangular mesh specified by the geometry g Point matrix p Edge matrix e and Triangle matrixt The triangular mesh is given by the mesh data p e andt Details on the mesh data representation can be found in the entry on i nit mesh pl e1 tl ul refinemesh g p e t u refines the mesh and also extends the function u to the new mesh by linear interpolation The number of rows in u should correspond to the number of columns in p and u1 has as many rows as there are points in p1 Each column of u is interpolated separately An extra input argument i t is interpreted as a list of subdomains to refine if it is a row vector or a list of triangles to refine if it is a column vector The default refinement method is regular refinement where all of the specified triangles are divided into four triangles of the same shape Longest edge
200. ostatics e Magnetostatics e AC Power Electromagnetics e Conductive Media DC e Heat Transfer e Diffusion These interfaces have dialog boxes where the PDE coefficients boundary conditions and solution are explained in terms of physical entities The application interfaces enable you to enter specific parameters such as Y oung s modulus in the structural mechanics problems Also visualization of the relevant physical variables is provided Several nontrivial examples are included in this manual Many examples are solved both by using the GUI and in command line mode The toolbox contains a number of demonstration M files They illustrate some ways in which you can write your own applications Can Extend the Functionality of the Toolbox The PDE Toolbox is written using MATLAB s open system philosophy There are no black box functions although some functions may not be easy to understand at first glance The data structures and formats are documented You can examine the existing functions and create your own as needed 1 Tutorial How Can I Solve 3 D Problems by 2 D Models The PDE Toolbox solves problems in two space dimensions and time whereas reality has three space dimensions The reduction to 2 D is possible when variations in the third space dimension taken to be z can be accounted for in the 2 D equation In some cases like the plane stress analysis the material parameters must be modified in the process of dime
201. ound Purpose Synopsis Description See Also 5 96 Write boundary condition specification file fid wbound bl mn fid wbound bl mn writes a Boundary M file with the name mn m The Boundary M file is equivalent to the Boundary Condition matrixb fid returns 1 if the file couldn t be written bl describes the boundary conditions of the PDE problem b is a Boundary Condition matrix See the entry on assemb for details The output file mn m is the name of a Boundary M file See the entry on pdebound decsg pdegeom pdebound wgeom wgeom Purpose Synopsis Description See Also Write geometry specification function fi d wgeom dl mn fid wgeom dl mn writes a Geometry M file with the name mn m The Geometry M file is equivalent to the Decomposed Geometry matrix dl fid returns 1 if the file couldn t be written dl isa Decomposed Geometry matrix The format of the Decomposed Geometry matrix is described in the entry on decsg The output file m m is the name of a Geometry M file The Geometry M file format is described in the entry on pdegeom decsg pdegeom wbound 5 97 A AC power electromagnetics 2 51 adaptive mesh refinement 1 22 4 26 adaptmesh 5 8 animation 2 18 3 31 3 36 Application command 3 11 application modes 2 35 Armijo Goldstein line search 4 22 Arnoldi algorithm 5 93 assema 5 13 assemb 5 14 assembling 4 7 5 18 assempde 5 18 Axes Limits
202. oundary conditions displays a dialog box in which you can specify the boundary condition for the selected boundary segments There are three different condition types e Generalized Neumann conditions where the boundary condition is determined by the coefficients q and g according to the following equation cVu qu 8 In the system cases q iS a 2 by 2 matrix and g is a 2 by 1 vector e Dirichlet conditions u is specified on the boundary The boundary condition equation is hu r where h is a weight factor that can be applied normally 1 In the system cases h is a 2 by 2 matrix andr is a 2 by 1 vector e Mixed boundary conditions system cases only which is a mix of Dirichlet and Neumann conditions q is a 2 by 2 matrix y is a 2 by 1 vector h isa 1 by 2 vector andr isa scalar PDE Toolbox M enus The following figure shows the boundary condition dialog box for the generic system PDE For boundary condition entries you can use the following variables in a valid MATLAB expression e The 2 D coordinates x and y e A boundary segment parameter s proportional to arc length s is 0 at the start of the boundary segment and increases to 1 along the boundary segment in the direction indicated by the arrow e Theoutward normal vector components nx andny If you need thetangential vector it can be expressed usingnx andny since tx ny and ty n e Thesolution u e Thetimet NOTE f the boundary condition is a
203. pplies to scalar PDEs If anongeneric application mode is selected application specific PDEs and parameters replace the standard PDE coefficients For a thorough description of the different application modes see the section Application Modes in the Examples chapter Each of the coefficientsc a f andd can be given as a valid MATLAB expression for computing coefficient values at the traingle centers of mass The following variables are available e x andy The x and y coordinates e u The solution e ux uy Thexand y derivatives of the solution et Thetime NOTE If the PDE coefficient is a function of the solution u or its derivatives ux and uy you must use the nonlinear solver If the PDE coefficient is a function of the timet you must choose a parabolic or hyperbolic PDE You can also enter the name of a user defined MATLAB function that accepts the arguments p t u ti me For an example type the function ci rclef 3 19 3 The Graphical User Interface 3 20 c can bea scalar or a 2 by 2 matrix The matrix c can be used to model e g problems with anisotropic material properties Ifc contains two rows they arethec and c elements of a 2 by 2 symmetric matrix If c contains three rows they are the c C1 2 and c elements of a 2 by 2 symmetric matrix C21 C 2 CHI e129 CS T ERI Ifc contains four rows they aretheG 1 G C 2 and G elements of the 2 by 2 matrix abo
204. presented by the column vector ele 0 E S CASA E Note that no lengths are stored for h or r Also for a scalar PDE the Dirichlet boundary condition u x2y is stored in the column vector Labels Ad ge SS O A A For a system N 2 with mixed boundary conditions M 1 Ay hipu r AED I T wl l s 4 122 82 assemb See Also the column looks like 2 1 q11 q21 q12 q22 gl 92 h11 h12 rl ql qi L q12 q22 Ol vas AVAE KLI aa h12 e A Wherel q11 1q21 denotelengths of the MATLAB text expressions andq11 q21 denote the actual expressions You can easily create your own examples by trying out pdet oo Enter boundary conditions by double clicking on boundaries in boundary mode and then export the Boundary Condition matrix tothe MATLAB main workspace by selecting the Export Decomposed Geometry Boundary Cond s option from the Boundary menu assempde pdebound 5 17 assempde Purpose Synopsis Description 5 18 Assemble the stiffness matrix and right hand side of a PDE problem K F B udJ assempde b p Ud assempde b p u assempde b p e t c a f u assempde b p e t c a f u0 u assempde b p e t c a f u0 ti me u assempde b p e t c a f time K F assempde b p e t c a f K F assempde b p e t c a f u0 K F assempde b p e t c a f Uu0 time K F assempde b p e t c a f u0 time sdl K F assempde b p e t c a f time K Fl assempde
205. ration and PDE solution Optional arguments are given as property name property value pairs The function produces a solution u to the elliptic scalar PDE problem V cVu au f on O or the elliptic system PDE problem V cC Vu au fon Q with the problem geometry and boundary conditions given by g and b The mesh is described by thep e andt The solution u is represented as the solution vector u Details on the representation of the solution vector can be found in the entry onassempde The algorithm works by solving a sequence of PDE problems using refined triangular meshes The first triangular mesh generation is obtained either as an optional argument toadapt mesh or by a call toi nit mesh without options The following generations of triangular meshes are obtained by solving the PDE problem computing an error estimate selecting a set of triangles based on the error estimate and then finally refining these triangles The solution to the PDE problem is then recomputed The loop continues until notriangles are selected by the triangle selection method or until the maximum number of trianglesis attained or until the maximum number of triangle generations has been generated g describes the decomposed geometry of the PDE problem g can either bea Decomposed Geometry matrix or the name of a Geometry M file The formats of the Decomposed Geometry matrix and Geometry M file are described in the entries on decsg and pdegeom respect
206. rations Each mesh point that is not located on an edge segment is moved toward the center of mass of the polygon formed by the adjacent triangles This process is repeated according to the setting of the0pt and Iter variables e When Opt is set toof f this process is repeated t er times default 1 e When Opt is set to mean the process is repeated until the mean triangle quality does not significantly increase or until the bound ter is reached default 20 e When Opt is set to mi n the process is repeated until the minimum triangle quality does not significantly increase or until the bound ter is reached default 20 Create a triangular mesh of the L shaped membrane first without jiggling and then jiggle the mesh p e t initmesh Ishapeg jiggle off q pdetriqlp t pdeplot p e t xydata q colorbar on xystyle flat pl jigglemesh p e t opt mean iter inf q pdetriq pl t pdeplot pl e t xydata g colorbar on xystyle flat initmesh pdetriq parabolic Purpose Synopsis Description Solve parabolic PDE problem ul parabolic u0 tlist b p e t c a f d ul parabolic u0 tlist b p e t c a f d rtol ul parabolic u0 tlist b p e t c a f d rtol atol ul parabolic u0 tlist K F B ud M ul parabolic u0 tlist K F B ud M rtol ul parabolic u0 tlist K F B ud M rtol atol ul parabolic u0 tlist g b p e t c a f d produces the solution tothe
207. raw mode by selecting Draw Mode fromthe Draw menu or by clicking on oneof the Draw mode icons to add another solid object Back in the Draw mode you are able to add change or delete solid objects and also to alter the set formula In addition to the recommended path of actions there are a number of shortcuts which allow you to skip over one or more steps In general the pdetool GUI adds the necessary steps automatically If you haven t yet defined a CSG model and leave the Draw mode with an empty model pdet oo creates an L shaped geometry with the default boundary condition and then proceeds to the action called for performing all the steps necessary If you arein the Draw mode and press the A button to initialize the mesh pdetool first decomposes the geometry using the current set formula and assigns the default boundary condition to the outer boundaries After that an initial mesh is created If you press the refine button to refine the mesh before the mesh has been initialized pdet oo first initializes the mesh and decomposes the geometry if you were still in the Draw mode If you press the button to solve the PDE and you have not yet created a mesh pdet oo initializes a mesh before solving the PDE If you select a plot type and choose to plot the solution pdet oo checks to see if thereis a solution tothe current PDE available If not pdet oo first solves the current PDE The solution is then displayed using the sele
208. rn off Toolbar Help Zoom Application Refresh Grid Grid Spacing Snap Axis Limits Axis E qual Turn off Toolbar Help Zoom Application Refresh Turn grid on off Adjust the grid spacing Turn the snap to grid feature on off Change the scaling of the drawing axes Turn the axis equal feature on off Turn off help texts for the toolbar buttons Turn zoom feature on off Select application mode Redisplay all graphical objects in the pdetool graphical user interface 3 The Graphical User Interface 3 10 Grid Spacing yj Grid Spacing X axis linear spacing JW Auto x axis extra ticks Y axis linear spacing JW Auto Y axis extra ticks Apply Done i In the Grid Spacing dialog box you can adjust the x axis and y axis grid spacing By default MATLAB s automatic linear grid spacing is used If you turn off the Auto checkbox the edit fields for linear spacing and extra ticks are enabled F or example the default linear spacing 1 5 0 5 1 5 can be changed to 1 5 0 2 1 5 1n addition you can add extra ticks so that the grid can be customized to aid in drawing the desired 2 D domain Extra tick entries can be separated using spaces commas semicolons or brackets Examples pi DES Ml 1 0 123 pi 4 Pressing the Apply push button applies the entered grid spacing pressing the Done button ends the Grid Spacing dialog PDE Too
209. s A column in bt corresponds to the column with the same index in gd A rowin bt corresponds to a minimal region index decsg dl bt dl1 bt1 msb decsg gd and dl bt dl 1 bt1 msb decsg gd sf ns return a second set of minimal regions dl 1 with a corresponding Boolean tablebt 1 This second set of minimal regions all have a connected boundary These minimal regions can be plotted by using MATLAB patch objects The second set of minimal regions have borders that may not have been induced by the original solid objects This occurs when two or more groups of solid objects have nonintersecting boundaries The calling sequences additionally return a sequence ms b of drawing commands for each second minimal region The first row contains the number of edge segment that bounds the minimal region The additional rows contain the sequence of edge segments from the Decomposed Geometry matrix that constitutes the bound If the index edge segment label is greater than the total number of edge segments it indicates that the total number of edge segments should be subtracted from the contents to get the edge segment label number and the drawing direction is opposite to the one given by the Decomposed Geometry matrix Geometry Description Matrix The Geometry Description matrix gd describes the CSG model that you draw using pdetool The current Geometry Description matrix can be made available to the MATLAB main workspace by selecting th
210. s calculated as follows Letcmax amax fmax andumax bethe maximum ofc a f andu respectively Let betheside of the smallest axis aligned square that contains the geometry Then scale max fmax 2 amax umax 2 cmax umax The scaling makes thet o parameter independent of the scaling of the equation and the geometry See Also adapt mesh pdej mps pdeadw orst Purpose Synopsis Description See Also Select triangles relative to the worst value bt pdeadworst p t c a f u errf wl evel bt pdeadworst p t c a f u errf wl evel returns indices of triangles to be refined in bt Used fromadapt mesh to select the triangles to be further refined The geometry of the PDE problem is given by the mesh datap andt Seethe entry oni nit mesh for more details c a andf arePDE coefficients See the entry on assempde for details u is the current solution given as a column vector See the entry on assempde for details errf isthe error indicator as calculated by pdej mps wl evel istheerror level relative to the worst error wl evel must be between O and 1 Triangles are selected using the criterion errf gt wl evel max errf adapt mesh pdej mps pdearcl Purpose Synopsis Description Examples See Also Interpolation between parametric representation and arc length pp pdearcl p xy s s0 s1 pp pdearcl p xy s s0 s1 returns parameter values for a parameterized curve corresponding to a given set of ar
211. s L and U sparseand the factorization numerically stable This factorization needs to be done only once in the beginning then x Ay is computed as x QU L PMy with one sparse matrix vector multiplication a permutation sparse forwards and backsubstitutions and a final renumbering 4 20 N onlinear Equations Nonlinear Equations The low level functions of the PDE Toolbox are aimed at solving linear equations Since many interesting computational problems are nonlinear the toolbox contains a nonlinear solver built on top of theassempde function The basic idea is to use Gauss N ewton iterations to solve the nonlinear equations Say you are trying to solve the equation r u V c u Vu a u u f u 0 In the FEM setting you solve the weak form of r u 0 Set as usual u x Y Uy multiply the equation by an arbitrary test function q integrate on the domain Q and use Green s formula and the boundary conditions to obtain 0 p U Elf cet U VO x VO x a x U 0 00 x dx J 20 UNO 009 0ds JU f fx UJO Cddx _ g a U O x ds Q dQ which has to hold for all indices i The residual vector p U can be easily computed as p U K M Q U F G where the matrices K M Q and the vectors F and G are produced by assembling the problem V c U Vu a U u f U 4 21 4 The Finite Element M ethod Assume that you have a guess U of the solution If U is close enough tothe exact
212. s N ewton method 1 21 1 Tutorial Small triangles are needed only in those parts of the computational domain where the error is large In many cases the errors are large in a small region and making all triangles small is a waste of computational effort Making small triangles only where needed is called adapting the mesh refinement to the solution An iterative adaptive strategy is the following For a given mesh form and solve the linear system KU F Then estimate the error and refine the triangles in which the error is large The iteration is controlled by adapt mesh and the error is estimated by pdej mps Although the basic equation is scalar systems of equations are also handled by the toolbox The interactive environment accepts u as a scalar or 2 vector function n command line mode systems of arbitrary size are accepted If c gt gt 0 and a gt 0 under rather general assumptions on the domain Q and the boundary conditions the solution u exists and is unique The FEM linear system has a unique solution which converges to u as the triangles become smaller The matrix K and theright hand side F make sense even when u does not exist or is not unique It is advisable that you devise checks to problems with questionable solutions 1 22 Using the Graphical User Interface Using the Graphical User Interface The PDE Toolbox Graphical User Interface The PDE Toolbox includes a complete graphical user interface GUI whi
213. s one row for each equation in the system There are two standard triangle selection methods in the toolbox pdeadworst andpdeadgsc pdeadworst selects triangles whereerrf exceedsa fraction default 0 5 of the the worst value andpdeadgsc selects triangles using a relative tolerance criterion Therefinement method is either ongest orregul ar Details on therefinement method can be found in the entry on ref i nemesh The adaptive algorithm can also solve nonlinear PDE problems For nonlinear PDE problems the Non in parameter must be set toon The nonlinear toleranceTo n nonlinear initial valueu0 nonlinear J acobian calculation J ac and nonlinear residual norm Nor m are passed to the nonlinear solver pdenonlin Seethe entry on pdenonlin for details on the nonlinear solver adaptmesh Examples Solve the Laplace equation over a circle sector with Dirichlet boundary conditions u cos 2 3atan2 y x along the arc and u 0 along the straight lines and compare to the exact solution We refine the triangles using the worst error criterion until we obtain a mesh with at least 500 triangles u p e t adaptmesh cirsg cirsb 1 0 0 maxt 500 tripick pdeadworst ngen inf x p 1 y p 2 exact x 2 y 7 2 1 3 cos 2 3 atan2 y x max abs u exact ans 0 0058 size t 2 ans 534 pdemesh p e t SEEKER TII ORAS SARPENADNS OSOS 5 11 adaptmesh The maximum absol
214. s over the unit square corresponding to its double second eigenvalue You may not have obtained the zero values along a diagonal of the square any linethrough the center of the square may have been computed This shows a general fact about multiple eigenvalues for symmetric matrices namely that any vector in the invariant subspace is equally valid as an eigenvector The two eigenfunctions Ug and Uy are orthogonal to each other if the dividing lines make right angles Check your solutions for that Actually the eigenvalues of the square can be computed exactly They are m24n2 n e g the double eigenvalue A and 219 is 10x which is pretty close to 100 If you compute the FEM approximation with only one refinement you would only find 16 eigenvalues and you obtain the wrong solution to the original problem You can of course check for this situation by computing the eigenvalues over a slightly larger range than the original problem 2 29 2 Examples You get some information from the printout in the MATLAB command window that is printed during the computation For this problem the algorithm computed a new set of eigenvalue approximations and tested for convergence every third step In the output you get the step number the time in seconds since the start of the eigenvalue computation and the number of converged eigenvalues with eigenvalues both inside and outside the interval counted Here is what MATLAB wrote
215. sary to satisfy the kinematic constraints described by the Dirichlet conditions The rest of this section details the treatment of the Dirichlet conditions and may be skipped on a first reading The PDE Toolbox supports two implementations of Dirichlet conditions The simplest is the Stiff Spring model so named for its interpretation in solid mechanics It has been mentioned in the section The Elliptic Equation on page 4 3 for the scalar case which is equivalent toa diagonal h matrix For the general case Dirichlet conditions hu r are approximated by adding a term L h hu h r to the equations KU F where L is a large number such as 10 times a representative size of the elements of K When this number is increased hu r will be more accurately satisfied but the potential ill conditioning of the modified equations will become more serious 4 11 4 The Finite Element M ethod The second method is also applicable to general mixed conditions with nondiagonal h and is free of the ill conditioning but is more involved computationally Assume that there are Np nodes in the triangulation Then the number of unknowns is N N Ny When Dirichlet boundary conditions fix some of the unknowns the linear system can be correspondingly reduced This is easily done by removing rows and columns when u values are given but here we must treat the case when some linear combinations of the components of u are given hu r Thes
216. sed as the set formula R1 R2 Leave the Draw mode and enter the Boundary mode by pressing the 0Q button and continue by selecting boundaries and specifying the boundary conditions U sing the Select All option from the E dit menu and then defining the Neumann condition eu 0 for all boundaries first is a good idea since that n leaves only the left most and right most boundaries to define individually The next step is to open the PDE Specification dialog box and enter the PDE coefficients 2 17 2 Examples The generic parabolic PDE that the PDE Toolbox solves is du V cV ds cCVu au f with initial values u0 u ty and the times at which to compute a solution specified in thearraytl ist For this case you have d 1 c 1 a O and f 0 Initialize the mesh by pressing the A button If you want you can refine the mesh by pressing the Refine button The initial values u0 0 and thelist of times is entered as the MATLAB array 0 0 5 5 They are entered into the Solve Parameters dialog box which is accessed by selecting Parameters from the Solve menu The problem can now be solved Pressing the button solves the heat equation at 11 different times from 0 to 5 seconds By default an interpolated plot of the solution i e the heat distribution at the end of the time span is displayed A more interesting way to visualize the dynamics of the heat distribution process is toanimatethe solution Tostart
217. solution an improved approximation UM 1 is obtained by solving the linearized problem an n 1 y _ n PEA ul Pu apo where a is a positive number It is well known that for sufficiently small pued lt pu and Py Pera ou is called a descent direction for p U where is thel norm The iteration is UOD U ap where a lt 1 is chosen as large as possible such that the step has a reasonable descent The Gauss N ewton method is local and convergence is assured only when UV is close enough to the solution In general the first guess may be outside the region of convergence To improve convergence from bad initial guesses a damping strategy is implemented for choosing a the Armijo Goldstein line search It chooses the largest damping coefficient out of the sequence 1 1 2 1 4 such that the following inequality holds PU p U ap 1 gt AUD which guarantees a reduction of the residual norm by at least 1 a1 2 Note that each step of the line search algorithm requires an evaluation of the residual p U ap An important point of this strategy is that when U approaches the solution then a gt 1 and thus the convergence rate increases If there is a solution to p U 0 the scheme ultimately recovers the quadratic convergence rate of the standard Newton teration 2 Itisnot necessary that p U 0 has a solution even if r u 0 has In this case the Gauss N
218. solution algorithms to work the mesh on the rectangle must be a regular mesh In this context it means that the first side of the rectangleis divided into N segments of length h the second into N segments of length hz and N 1 by N 1 points are introduced on the regular grid thus defined The triangles are all congruent with sides h h and a right angle in between The Dirichlet boundary conditions are eliminated in the usual way and the resulting problem for the interior nodes is Kv F If the interior nodes are numbered from left to right and then from bottom to top the K matrix is block tridiagonal The N 1 diagonal blocks here called T are themselves tridiagonal N 1 by N 1 matrices with 2 h h h gt h on the diagonal and h h on the subdiagonals The subdiagonal blocks here called are h h times the unit N 1 matrix The key to the solution of the problem Kv F is that the problem Tw f is possible to solve using the discrete sinetransform Let S be the N 1 by N 1 matrix with Sij sin rij N Then SITS A where A is a diagonal matrix with diagonal entries 2 h h gt h hy 2h h cos ri N 1 w SA 1S7 f but multiplying with S is nothing morethan takingthe discrete sinetransform and multiplying with S is the same as taking the inverse discrete sine transform The discrete sine transform can be efficiently calculated using the fast Fourier transform on a sequence of length 2N
219. splay the node number on the information line Press the middle mouse button or use the left mouse button and the Shift key to display the triangle number on the information line In Plot mode you can use the mouse to display the numerical value of the plotted property at the position where you click Press the left mouse button to display the triangle number and the value of the plotted property on the information line The information remains on the information line until you release the mouse button 1 35 1 Tutorial Entering Parameter Values as MATLAB Expressions When entering parameter values e g as a function of x and y the entered string must bea MATLAB expression to be evaluated for x and y defined on the current mesh e xand y are MATLAB row vectors For example the function 4 xy should be entered as 4 x y and not as 4 x y which normally is not a valid MATLAB expression Using PDE Toolbox version 1 0 Model M files You can convert M odd M files created using version 1 0 of the PDE Toolbox for MATLAB 4 2 for use with PDE Toolbox version 1 0 2 and MATLAB 5 Theold Model M files cannot be used directly in PDE Toolbox 1 0 2 To convert your old Model M files use the conversion utility pdemdl cv To convert a Model M file called model 42 m toa PDE Toolbox 1 0 2 compatible Model M filecalled mode 5 m type the following at the MATLAB command line pdemdlcv model 42 model 5 1 36 Using Command line Functio
220. t important functions They do not offer any additional functionality all toolbar buttons functions are also available using menu items The toolbar consists of three different parts the five left most buttons for draw mode functions the next six buttons for different boundary mesh solution and plot functions and finally the right most button for activating the zoom feature Of a oje gt 30 pol A aj A The draw mode buttons represent from left to right Draw a rectangle square starting at a corner Using the left mouse button click and drag to create a rectangle Using the right mouse button or Control click click and drag to create a square Draw a rectangle square starting at the center Using the left mouse button click and drag to create a rectangle Using the right mouse button or Control click click and drag to create a square Draw an ellipse circle starting at the perimeter Using the left mouse button click and drag to create an ellipse Using the right mouse button or Control click click and drag to create a circle Draw an ellipse circle starting at the center Using the left mouse button click and drag to create an ellipse Using the right mouse button or Control click click and drag to create a circle Draw a polygon Click and drag to create polygon edges Y ou can close the polygon by pressing the right mouse button Clicking at the starting vertex also closes the polygon The draw mode butto
221. t of pdenonl i n the coefficients can depend on u The coefficients cannot depend ont the time A complete listing of all format options can be found in the entry onassempde pdenonlin The solver can be fine tuned by setting some of the options described below Property Name Property Value Default Description Jacobian fixed lumped full fixed Approximation of J acobian U0 string or numeric 0 Initial solution guess Tol positive scalar le 4 Residual size at termination Maxlter positive integer 25 Maximum Gauss N ewton iterations MinStep positive scalar 1 2 16 Minimum damping of search direction Report on off of f Print convergence information Norm string or numeric Inf Residual norm There are three methods currently implemented to compute the J acobian e Numerical evaluation of the full J acobian based on the sparse version of the function numjac e A lumped approximation described in Chapter 4 The Finite Element Method based on the numerical differentiation of the coefficients e A fixed point iteration matrix where the J acobian is approximated by the stiffness matrix Select the desired method by setting the acobi an property tof ul lumped or fi xed bearing in mind that the more precise methods are computationally more expensive U0 is the starting guess that can be given as an expression a generic scalar or a vector By default it is set to 0 but this is useless in problems such as V 1 uV u
222. tax of each of the text expressions must be according to the above item The number of expressions in the sequence must equal the number of subdomain labels in the triangle list t This number can be checked by typing max t 4 The name of a user defined MATLAB function that accepts the arguments p t u ti me u andti me areempty matrices if the corresponding parameter isnot passed toassempde p andt are mesh data u istheu0 input argument andt isthet i me input argument toassempde IftimeisNaN andthe function depends on time the function must return a matrix of correct size containing Na Ns in all positions assempde We refer to the matrices above as Coeffi cient matrix and the user defined MATLAB function as Coefficient M file Ifc contains two rows with data according to any of the above tems they are the cy 1 and cp 2 elements of the 2 by 2 diagonal matrix Ifc contains four rows they arethec j Cp 1 2 and Cp 2 elements of a 2 by 2 matrix PDE Coefficients for System Case Let N bethe dimension of the PDE system Now c is an N by N by 2 by 2 tensor a an N by N matrix and f a column vector of length N The elements Gijki aij Ajj and f of a d and f are stored row wisein the MATLAB matrices c a d andf Each row in these matrices is similar in syntax tothe scalar case There is one difference however At the point of evaluation of MATLAB text expressions the variablesu ux and uy contain matrices
223. ted as a combination u x ee Udi where p are some special piecewise linear basis functions and U are scalar coefficients Choose q like a tent such that it has the height 1 at the nodei and the height 0 at all other nodes For any fixed v the FEM formulation yields an algebraic equation in the unknowns U You want to determine N unknowns so you need N different instances of v What better candidates than v 9 j 1 2 N You finda linear system KU F where the matrix K and theright hand side F contain integrals in terms of the test functions 9 and the coefficients defining the problem c a f q and g The solution vector U contains the expansion coefficients of u which are also the values of u at each node x since u x U If the exact solution u is smooth then FEM computes u with an error of the same size as that of the linear interpolation It is possible to estimate the error on each triangle using only u and the PDE coefficients but not the exact solution u which in general is unknown The PDE Toolbox provides functions that assemble K and F This is done automatically in the graphical user interface but you also have direct access to the FEM matrices from the command line function assempde To summarize the FEM approach is to approximate the PDE solution u bya piecewise linear function 4 4 is expanded in a basis of test functions 9 and the residual is tested against all the basis functions
224. ten useful to decompose it into the union of more subdomains of simpler structure Such structures are often introduced by pdetool Assume now that Q is the disjoint union of some subdomains Q4 Q5 Q Then you could renumber the nodes of a mesh on Q such that the indices of the nodes of each subdomain are grouped together while all the indices of nodes common to two or more subdomains come last Since K has nonzero entries only at the lines and columns that are indices of neighboring nodes the stiffness matrix is partitioned as follows K 0 0 BI T 0 K OB 0 0 KB n B B B C n T n while the right hand side is The PDE Toolbox routineassempde can assemble the matrices Kj Bj fj and C separately You have full control over the storage and further processing of these matrices Furthermore the structure of the linear system Ku F is simplified by decomposing K into the partitioned matrix above Examples of Elliptic Problems Now consider the geometry of the L shaped membrane Y ou can plot the geometry of the membrane by typing pdegplot Ishapeg Noticethe borders between thesubdomains There arethree subdomains Thus the matrix formulas with n 3 from above can be used Now generate a mesh for the geometry p e t initmesh Ishapeg p e t refinemesh shapeg p e t p e t refinemesh shapeg p e t So for this case with n 3 you have T K 0 O Bi uy fi 0 K
225. the dialog box you can enter any device option that is available for MATLAB s pri nt command The default device option is dps PostScript for black and white printers The paper orientation can beset toportrait landscape ortal l and you can print toa printer or tofile PDE Toolbox M enus Edit Menu Edit Options Draw Undo lt Ctrl gt z Cut lt Ctrl gt x Copy lt Ctrl gt c Paste lt Ctrl gt v Clear lt Ctrl gt r Select All lt Ctrl a Undo Undo the last line when drawing a polygon Cut Move the selected solid objects to the Clipboard Copy Copy the selected objects to the Clipboard leaving them intact in their original location Paste Copy the contents of the Clipboard to the current CSG model Clear Delete the selected objects Select All Select all solid objects in the current CSG model Also select all outer boundaries or select all subdomains 3 The Graphical User Interface Paste displays a dialog box for pasting the contents of the Clipboard on to the current CSG model The Clipboard contents can be repeatedly pasted adding a specified x and y axis displacement to the positions of the Clipboard objects Using the default values zero displacement and one repetition the Clipboard contents is inserted at its original position PDE Toolbox M enus Options Menu Options Draw Boundar x Grid Grid Spacing x Snap Axes Limits Axes Equal Tu
226. the main workspace Then makea script M fileor type the following commands at the MATLAB prompt gt hmnewplot hf get h Parent set hf Renderer zbuffer axis tight set gca DataAspectRatio 1 1 1 axis off M moviein 10 hf gt maxu max abs u col ormap cool for j 1 10 ur real exp j 2 pi 10 sqrt 1 u pdeplot p e t xydata ur colorbar off mesh off caxis maxu maxu axis tight set gca DataAspectRatio 1 1 1 axis off M j getframe end movie hf M 50 pdedemo2 contains a full command line demonstration of the scattering problem 2 9 2 Examples A Minimal Surface Problem In many problems the coefficients c a and f donot only depend on x and y but also on the solution u itself Consider the equation V 0 1 Vu on the unit disk Q y x y lt 1 with u x on 20 This problem is nonlinear and cannot be solved with the regular elliptic solver Instead the nonlinear solver pdenonl in is used Let s solve this minimal surface problem using thepdetoo GUI and command line functions Using the Graphical User Interface Make sure that the application mode in thepdetoo GUI is set to Generic Scalar The problem domain is simply a unit circle Draw it and move to the Boundary mode to define the boundary conditions Use Select All from the E dit menu to select all boundaries Then double click on a boundary to open the Boun
227. the type of object that you want to use by clicking on the corresponding button or by using the Draw pull down menu Combine the solid objects and the set algebra to build the desired CSG model 1 31 1 Tutorial 2 Savethe geometry to a model file The model file is an M file so if you want to continue working using the same geometry at your next PDE Toolbox session simply type the name of the model file at the MATLAB prompt The pdetool GUI then starts with the model file s solid geometry loaded If you save the PDE problem at a later stage of the solution process the model file also contains commands to recreate the boundary conditions the PDE coefficients and the mesh 3 Movetothenext step in the PDE solving process by pressing the dQ button The outer boundaries of the decomposed geometry are displayed with the default boundary condition indicated If the outer boundaries do not match the geometry of your problem re enter the Draw mode You can then correct your CSG model by adding removing or altering any of the solid objects or change the set formula used to evaluate the CSG model NOTE The set formula can only be edited while you are in the Draw mode If the drawing process resulted in any unwanted subdomain borders remove them by using the Remove Subdomain Border or Remove All Subdomain Borders option from the Boundary menu Y ou can now define your problem s boundary conditions by selecting the boundary to chan
228. tional steps in the process of modeling and solving your PDE are then saved to the same M file This concludes the drawing part You can now define the boundary conditions for the outer boundaries Enter the Boundary mode by pressing the dQ icon or by selecting Boundary Mode from the Boundary menu You can now remove subdomain borders and define the boundary conditions The gray edge segments are subdomain borders induced by the intersections of the original solid objects Borders that do not represent borders between e g areas with differing material properties can be removed From the Boundary menu select the Remove All Subdomain Borders option All borders are then removed from the decomposed geometry The boundaries are indicated by colored lines with arrows The color reflects the type of boundary condition and the arrow points towards the end of the boundary segment The direction information is provided for the case when the boundary condition is parameterized along the boundary The boundary condition can also bea function of x and y or simply a constant By default the boundary condition is of Dirichlet type u 0 on the boundary Dirichlet boundary conditions are indicated by red color The boundary conditions can also be of a generalized Neumann blue or mixed green type For scalar u however all boundary conditions are either of Dirichlet or the generalized Neumann type You select the boundary conditions that you want to c
229. tions and applications of FEM to problems in structural mechanics see 1 Example Consider a steel platethat is clamped along a right angleinset at the lower left corner and pulled along a rounded cut at the upper right corner All other sides are free The steel plate has the following properties Dimension 1 by 1 meters thickness 1mm inset is 1 3 by 1 3 meters The rounded cut runs from 2 3 1 to 1 2 3 Young s modulus 196 10 MN m Poisson s ratio 0 31 The curved boundary is subjected to an outward normal load of 500 N m We need to specify a surface traction we therefore divide by the thickness 1 mm thus thesurfacetractions should beset to0 5 MN m We will usethe force unit MN in this example We want to compute a number of interesting quantities such as the x and y direction strains and stresses the shear stress and the von Mises effective stress Using the Graphical User Interface Using the pdet 00 GUI the first thing to do is to select the application mode For this problem use Structural Mechanics Plane Stress The CSG model can be made very quickly by drawing a polygon with corners in x 0 2 3 1 1 1 3 1 3 0 andy 1 1 2 3 0 0 1 3 1 3 anda circle with center inx 2 3 y 2 3 and radius1 3 2 39 2 Examples The polygon is normally labeled P1 and the circle C1 and the CSG model of the steel plate is simply P14C1 Next select Boundary Mode to specify the boundary conditions First remove a
230. tor of values on the nodes of the current mesh u t0 The initial value u t for the hyperbolic PDE problem Y ou can use the same formats as for u tO Relative tolerance Relative tolerance parameter for the ODE solver that is used for solving the time dependent part of the hyperbolic PDE problem Absolute tolerance Absolute tolerance parameter for the ODE solver that is used for solving the time dependent part of the hyperbolic PDE problem PDE Toolbox M enus Solve Parameters Eigenvalue search range o 100 OK Cancel Solve Parameters dialog box for eigenvalue PDEs all e Eigenvalue problems For the eigenvalue PDE the only solve parameter is the Eigenvalue search range a two element vector defining an interval on the real axis as a search range for the eigenvalues The left hand side can be Inf Examples 0 100 Inf 50 3 29 3 The Graphical User Interface 3 30 Plot Menu Plot Solution Display a plot of the solution Parameters Open dialog box for plot selection Export Movie If a movie has been recorded the movie matrix M is exported to the main workspace Parameters The Plot Selection dialog box Parameters opens a dialog box containing options controlling the plotting and visualization PDE Toolbox M enus The upper part of the dialog box contains four columns Plot type far left contains a row of six different plot types which can be
231. trix or the name of a Boundary M file The formats of the Boundary Condition matrix and Boundary M file are described in the entries on assemb andpdebound respectively The geometry of the PDE problem is given by the mesh data p e andt Details on the mesh data representation can be found in the entry oni nit mesh The optional list of subdomain labels s dl restricts the assembly process to the subdomains denoted by the labels in the list The optional input arguments u 0 andti me are used for the nonlinear solver and time stepping algorithms respectively The tentative input solution vector u0 has the same format as u PDE Coefficients for Scalar Case The coefficients c a and fin the scalar PDE problem can be represented in the MATLAB variablesc a andf in the following ways e A constant e A row vector of values at the triangle centers of mass e A MATLAB text expression for computing coefficient values at the triangle centers of mass The expression is evaluated in a context wherethe variables x y 5d u ux uy andt are row vectors representing values at the triangle centers of mass t is a scalar The row vectors contain x and y coordinates subdomain label solution x and y derivatives ofthesolution andtime u ux and uy can only be used if u0 has been passed toa s s e mp de The same applies tothe scalar t which is passed toassempde asti me A sequence of MATLAB text expressions separated by exclamation marks The syn
232. trostatics the electrostatic scalar potential V is related tothe electric field E by E VV and using one of Maxwdl s equations V D p and the relationship D E we arrive at the Poisson equation _V eVV p where e is the coefficient of die ectricity and p is the space charge density NOTE should really be written as ego where e is the coefficient of dielectricity or permittivity of vacuum 8 854 10 12 farad meter and e is the relative coefficient of dielectricity that varies among different dielectrics 1 00059 in air 2 24 in transformer oil etc Using the PDE Toolbox application mode Electrostatics you can solve electrostatic problems modeled by the equation above The PDE Specification dialog box contains entries for e and p The boundary conditions for electrostatic problems can be of Dirichlet or Neumann type For Dirichlet conditions the electrostatic potential V is specified on the boundary For Neumann conditions the surface charge n eVV is specified on the boundary For visualization of the solution to an electrostatic problem the plot selections include the electrostatic potential V the electric field E and the electric displacement field D For a more in depth discussion of problems in electrostatics see 2 2 43 2 Examples Example Let s consider the problem of determining the electrostatic potential in an air filled quadratic frame bounded by a square with side length
233. tton or the PDE Specification menu item from the PDE menu NOTE This step can be performed at any time prior to solving the PDE since it is independent of the CSG model and the boundaries If the PDE coefficients are material dependent they are entered in the PDE mode by double clicking on the different subdomains 7 Solvethe PDE by pressing the button or by selecting Solve PDE from the Solve menu If you don t want an automatic plot of the solution or if you want to change the way the solution is presented you can do that from the Plot Selection dialog box prior to solving the PDE You open the Plot Selection dialog box by pressing the button with the 3 D solution plot icon or by selecting the Parameters menu item from the Plot menu 8 Now from here you can choose one of several alternatives Export the solution and or the mesh to MATLAB s main workspace for further analysis Visualize other properties of the solution Change the PDE and recompute the solution Change the mesh and recompute the solution If you select Initialize Mesh the mesh is initialized if you select Refine Mesh the current mesh is refined From the Mesh menu you can also jiggle the mesh and undo previous mesh changes Change the boundary conditions To return to the mode where you can se lect boundaries usethe dQ button or the Boundary Mode option from the Boundary menu 1 Tutorial e Change the CSG model You can re enter the d
234. ur disposal to solve the linear system and further process the solution More elaborate applications make use of the FEM specific information returned by the different functions of the toolbox Therefore we quickly summarize the theory and technique of FEM solvers to enable advanced applications to make full use of the computed quantities FEM can besummarized in the following sentence Project the weak form of the differential equation onto a finitedimesional function space The rest of this section deals with explaining the above statement Westart with the weak form of the differential equation Without restricting the generality we assume generalized Neumann conditions on the whole boundary since Dirichlet conditions can be approximated by generalized Neumann conditions In the simple case of a unit matrix h setting g qr and then letting q gt yields the Dirichlet condition because division with a very large q cancels the normal derivative terms The actual implementation is different since the above procedure may create conditioning problems The mixed boundary condition of the system case requires a more complicated treatment detailed in the section The Elliptic System Assume that u is a solution of the differential equation Multiply the equation with an arbitrary test function v and integrate on Q V cVu v auv dx fv dx Q Q Integrate by parts i e use Green s formula to obtain The boundary integra
235. used for visualization Color Visualization of a scalar property using colored surface objects Contour Visualization of a scalar property using colored contour lines The contour lines can also enhance the color visualization when both plot types Color and Contour are checked The contour lines are then drawn in black Arrows Visualization of a vector property using arrows Deformed mesh Visualization of a vector property by deforming the mesh using the vector property The deformation is automatically scaled to 10 of the problem domain This plot type is primarily intended for visualizing x and y displacements u and v for problems in structural mechanics If no other plot type is selected the deformed triangular mesh is displayed Height 3 D plot Visualization of a scalar property using height Z axis in a 3 D plot 3 D plots are plotted in separate figure windows If the Color and Contour plot types are not used the 3 D plot is simply a mesh plot Y ou can visualize another scalar property simultaneously using Color and or Contour which results in a 3 D surface or contour plot Animation Animation of time dependent solutions to parabolic and hyperbolic problems If you check this option the solution is recorded and then animated in a separate figure window using MATLAB s movie function 3 31 3 The Graphical User Interface 3 32 A colorbar is added to the plots to map the colors in the plot to the magnitude
236. using a contour plot and by plotting the heat flux vector field using arrows ial PDE Toolbox Untitled File Edit Options Draw Boundary PDE Mesh Solve Plot Window Help Time 0 1 Color T Vector field q T 35 A 1 r T T T T 3H 4 0 35 25 al 0 3 eb gt B 4 0 25 1 5 8 4 0 2 ik 4 0 15 asp 4 0 1 ob 4 0 05 05 i i i i i i i i i i 0 1 0 05 1 15 2 25 3 35 4 Info Select a different plot or change mode to alter PDE mesh or boundary conditions Exit Visualization of the temperature and the heat flux 2 60 Application Modes Diffusion Since heat transfer is a diffusion process the generic diffusion equation has the same structure as the heat equation dc 2 _V DVc En DVc Q where c is the concentration D is the diffusion coefficient and Q is a volume source The diffusion process may be anisotropic in which case D is a 2 by 2 matrix The boundary conditions can be of Dirichlet type where the concentration on the boundary is specified or of Neumann type where the flux n DV c is specified It is also possible to specify a generalized Neumann condition It is defined by n DV c qc g where q is a transfer coefficient Visualization of the concentration its gradient and the flux is available from the Plot Selection dialog box 2 61 2 Examples 2 62 The Graphical User Interface PDE Toolbox Menus
237. ute error is 0 0058 with 534 triangles We test how many refinements we have to use with an uniform triangle net p e t initmesh cirsg p e t refinemesh cirsg p e t usassempde cirsb p e t 1 0 0 x p 1 y p 2 exact x 2 y 2 1 3 cos 2 3 atan2 y x max abs u exact ans 0 0085 size t 2 ans 1640 p e t refinemesh cirsg p e t us assempde cirsb p e t 1 0 0 x p 1 y p 2 exact x 2 y 2 1 3 cos 2 3 atan2 y x max abs u exact ans 0 0054 size t 2 ans 6560 pdemesh p e t Thus with uniform refinement we need 6560 triangles to achieve better absolute error than what we achieved with the adaptive method Note that the error is reduced only by 0 6 when the number of elements in quadrupled by the uniform refinement For a problem with regular solution we expect a O h error but this solution is singular since u r at the origin Diagnostics Upon termination one of the following messages is displayed e Adaption completed This means that theTri pi ck function returned zero triangles to refine e Maximum number of triangles obtained e Maximum number of refinement passes obtained See Also initmesh refinemesh assempde pdeadgsc pdeadworst pdej mps 5 12 assema Purpose Synopsis Description See Also Assemble area integral contributions assem assem assem assem assem
238. ution specified in the arrayt ist For this case you haved 1 c 1 a 0 and f 0 The initial value u0 0 and the list of times tl ist is set tothe MATLAB array0 0 5 5 2 19 2 Examples To compute the solution call parabolic u parabolic 0 0 0 5 5 crackb p e t 1 0 0 1 The solution u created this way a matrix with 11 columns where each column corresponds to the solution at the 11 points in time 0 0 5 4 5 5 0 Let s plot the solution at t 5 0 seconds using interpolated shading anda hidden mesh Use the hot colormap pdeplot p e t xydata u 11 mesh off colormap hot 2 20 Examples of Parabolic Problems Heat Distribution in Radioactive Rod This heat distribution problem is an example of a 3 D parabolic PDE problem that is reduced to a 2 D problem by using cylindrical coordinates Consider a cylindrical radioactive rod At the left end heat is continuously added The right end is kept at a constant temperature At the outer boundary heat is exchanged with the surroundings by transfer At the same time heat is uniformly produced in the whole rod dueto radioactive processes Assume that the initial temperature is zero This leads to the following problem du pes V kVu f where p is the density C is the rod s thermal capacity k is the thermal conductivity and f is the radioactive heat source The density for this metal rod is 7800 kg m the thermal capacity is 500 Ws
239. ve The available types of PDEs are e Elliptic The basic form of the elliptic PDE is V cVu au f The parameter d does not apply to the elliptic PDE Parabolic The basic form of the parabolic PDE is Y cVu au f with initial values u0 Uu to Hyperbolic The basic form of the hyperbolic PDE is ia f ot with initial values ud u ty and ut 0 ty Eigenmodes The basic form of the PDE eigenvalue problem is V cVu au Adu The parameter f does not apply to the eigenvalue PDE PDE Toolbox M enus In the system case c is a rank four tensor which can be represented by four 2 by 2 matrices c11 12 21 andc22 They can beentered as one two three or four rows see the scalar case above a and d are 2 by 2 matrices and f is a 2 by 1 vector The PDE Specification dialog box for the system case is shown below 3 21 3 The Graphical User Interface Mesh Menu Mesh Solve Plot Window Help Mesh Mode Initialize Mesh lt Ctrl gt i Refine Mesh lt Ctrl m Jiggle Mesh Undo Mesh Change Display Triangle Quality Show Node Labels Show Triangle Labels Parameters Export Mesh Mesh Mode Initialize Mesh Refine Mesh J iggle Mesh Undo Mesh Change Display Triangle Quality 3 22 Enter mesh mode Build and display an initial triangular mesh Uniformly refine the current triangular mesh J iggle the mesh Undo the last mesh change All mesh generations are saved
240. ve definite The generalized eigenvalue problem KU AMU is now solved by the Arnoldi algorithm applied to a shifted and inverted matrix with restarts until all eigenvalues in the user specified interval have been found Let us describe how this is donein more detail You may want to look at the example provided in the section Examples of Eigenvalue Problems on page 2 27 where an actual run is reported First a shift u is determined close to where we want to find the eigenvalues When both K and M are positive definite it is natural totakeu 0 and get the smallest eigenvalues in other cases take any point in theinterval Ib ub where eigenvalues are sought Subtract uM from the eigenvalue equation and get K uM U A WMU Then multiply with the inverse of this shifted matrix 4 17 4 The Finite Element M ethod and get eri K uM MU This is a standard eigenvalue problem AU 6U with the matrix A K uM M and eigenvalues 9 Lu 1 1 n Thelargest eigenvalues 6 of the transformed matrix A now correspond tothe eigenvalues A u 1 6 of the original pencil K M closest to the shift u TheArnoldi algorithm computes an orthonormal basis V where the shifted and inverted operator A is represented by a Hessenberg matrix H AV VjH E The subscripts mean that Vj and Ej havej columns and Hj j hasj rows and columns When no subscripts are used we deal with vectors and matrices of size n Some
241. ving particular problems involving different kinds of PDEs geometries and boundary conditions and covering a range of different applications 1 17 1 Tutorial 1 18 Basics of The Finite Element Method The solutions of simple PDEs on complicated geometries can rarely be expressed in terms of elementary functions You are confronted with two problems First you need to describe a complicated geometry and generate a mesh on it Then you need to discretize your PDE on the mesh and build an equation for the discrete approximation of the solution Thepdet ool graphical user interface provides you with easy to use graphical tools to describe complicated domains and generatetriangular meshes It alsodiscretizes PDEs finds discrete solutions and plots results You can access the mesh structures and the discretization functions directly from the command line or M file and incorporate them into specialized applications Below is an overview of the Finite Element Method FEM The purpose of this presentation is to get you acquainted with the elementary FEM notions Here you find the precise equations that are solved and the nature of the discrete solution Different extensions of the basic equation implemented in the PDE Toolbox are presented A more detailed description can be found in Chapter 4 The Finite Element Method You start by approximating the computational domain Q with a union of simple geometric objects in this case triangles
242. with respect to u K and M designate stiffness and mass matrices and their indices designate the coefficients with respect to which they are assembled At each Gauss Newton iteration the nonlinear solver assembles the matrices corresponding tothe equations 4 24 N onlinear Equations V cvu a f u 0 and V cVu a u 0 and then produces the approximate J acobian The differentiations of the coefficients are done numerically In the general setting of elliptic systems the boundary conditions are appended to the stiffness matrix to form the full linear system GEA where the coefficients of K and F may depend on the solution The lumped approach approximates the derivative mapping of the residual by fro H 0 The nonlinearities of the boundary conditions and the dependencies of the coefficients on the derivatives of U are not properly linearized by this scheme When such nonlinearities are strong the scheme reduces to the fix point iteration and may converge slowly or not at all When the boundary conditions are linear they do not affect the convergence properties of the iteration schemes In the Neumann case they are invisible H is an empty matrix and in the Dirichlet case they merely state that the residual is zero on the corresponding boundary points 4 25 4 The Finite Element M ethod 4 26 Adaptive Mesh Refinement Thetoolbox has a function for global uniform mesh refinement It divides ea
243. xamples See Also Shorthand command for surface plot pdesurf p t u pdesurf p t u plots a 3 D surface of PDE node or triangle data Ifu isa column vector node data is assumed and continuous style and interpolated shading are used If u is a row vector triangle data is assumed and discontinuous style and flat shading are used h pdesurf p t u additionally returns handles to the drawn axes objects For node data this command is just shorthand for the call pdeplot p t xydata u xystyle interp zdata u zstyle continuous colorbar off and for triangle data it is pdeplot p t xydata u xystyle flat zdata u zstyle discontinuous colorbar off If you want to have more control over your surface plot usepdep ot instead of pdesurf Surface plot of the solution to the equation Au 1 over the geometry defined by the L shaped membrane Use Dirichlet boundary conditions u 0 on Q p e t initmesh Ishapeg p e t refinemesh Ishapeg p e t Uu assempde Ishapeb p e t 1 0 1 pdesurf p t u pdeplot pdecont pdemesh 5 77 pdetool Purpose Synopsis Description 5 78 PDE Toolbox graphical user interface GUI pdetool pdetool action flag pdetool provides the graphical user interface GUI for the PDE Toolbox Call pdetool without arguments to start the application Y ou should not call pdetoo with arguments The GUI helps you to draw the 2 D d
244. ystem 2 MoS KU F dt with the initial conditions U 0 ug x UNO vo xj Vi after we eliminatethe unknowns fixed by Dirichlet boundary conditions As before the stiffness matrix K and the mass matrix M are assembled with the aid of the function assempde from the problems V cVu au f and V 0Vu du 0 4 16 The Eigenvalue Equation The Eigenvalue Equation The basic eigenvalue problem handled by the PDE Toolbox is V cVu au Adu where a is an unknown complex number In solid mechanics this is a problem associated with wave phenomena describing e g the natural modes of a vibrating membrane In quantum mechanics A is the energy level of a bound state in the potential well a x The numerical solution is found by discretizing the equation and solving the resulting algebraic eigenvalue problem Let us first consider the discretization Expand u in the FEM basis multiply with a basis element and integrate on the domain Q This yields the generalized eigenvalue equation KU AMU where the mass matrix corresponds to the right hand side i e M j dOd The matrices K and M are produced by calling ass e ma for the equations V cVu au O and V OVu du 0 Notethat in the most common case when the function d x is positive the mass matrix M is positive definite symmetric Likewise when c x is positive and we have Dirichlet boundary conditions the stiffness matrix K is also positi

Partial Differential Equation Toolbox User's Guide

Contents

Download Pdf Manuals

Related Search

Related Contents