MAKE A MEME
View Large Image
| View Original: | Finite element triangulation.svg (652x652) | |||
| Download: | Original | Medium | Small | Thumb |
| Courtesy of: | commons.wikimedia.org | More Like This | ||
| Keywords: Finite element triangulation.svg Illustration of the en Finite element method the triangulation of the domain self-made with en Matlab 02 18 15 June 2007 UTC Oleg Alexandrov Source code MATLAB The triangulation used in this code and shown above is created with the http //www cs cmu edu/~quake/triangle html Triangle mesh generator <source lang matlab > Solve the problem -\Delta u + c u f with Dirichlet boundary conditions The domain is the disk centered at the origin of radius 1 Its triangulation is read from the end of this file function main c 0; a parameter in the equation see above white 0 99 1 1 1; blue 0 129 205/256; green 0 200 70/256; load the triangulation from the end of the file dummy_arg 0; P get_points dummy_arg ; T get_triangles dummy_arg ; find the number of points and the number of triangles Np k size P ; Nt k size T ; Nb 30; first Nb points in P are on the boundary plot the triangulation lw 1 4; figure 1 ; clf; hold on; axis equal; axis off; for l 1 Nt i T l 1 ; j T l 2 ; k T l 3 ; plot the three edges in each triangle plot P i 1 P j 1 P i 2 P j 2 'linewidth' lw 'color' blue plot P i 1 P k 1 P i 2 P k 2 'linewidth' lw 'color' blue plot P j 1 P k 1 P j 2 P k 2 'linewidth' lw 'color' blue end a hack to deal with bounding box issues s 1 05; plot -s -s ' ' 'color' white ; plot s s ' ' 'color' white ; save as eps and svg needs the plot2svg function saveas gcf 'triangulation eps' 'psc2' plot2svg 'Finite_element_triangulation svg' ; right-hand side f inline '5-x 2-y 2' 'x' 'y' ; f inline '4+0 x 2+0 y 2' 'x' 'y' ; the values of f at the nodes in the triangulation F f P 1 P 2 ; RHS 0 F; RHS RHS Nb+1 Np ; will solve A U RHS an empty sparse matrix of size Np by Np A sparse Np-Nb Np-Nb ; iterate through triangles for l 1 Nt fill in the matrix i T l 1 ; j T l 2 ; k T l 3 ; pii pij pik pjj pjk pkk gii gij gik gjj gjk gkk calc_elems i j k P i P j P k ; One has to consider the cases when a given vertex is on the boundary or not If yes it can't be included in the matrix if i > Nb A i-Nb i-Nb A i-Nb i-Nb + gii + c pii; end if i > Nb j > Nb A i-Nb j-Nb A i-Nb j-Nb + gij + c pij; A j-Nb i-Nb A i-Nb j-Nb ; end if i > Nb k > Nb A i-Nb k-Nb A i-Nb k-Nb + gik + c pik; A k-Nb i-Nb A i-Nb k-Nb ; end if j > Nb A j-Nb j-Nb A j-Nb j-Nb + gjj + c pjj; end if j > Nb k > Nb A j-Nb k-Nb A j-Nb k-Nb + gjk + c pjk; A k-Nb j-Nb A j-Nb k-Nb ; end if k > Nb A k-Nb k-Nb A k-Nb k-Nb + gkk + c pkk; end add the appropriate contributions to the right-hand terms if i > Nb RHS i-Nb RHS i-Nb + F i pii + F j pij + F k pik; end if j > Nb RHS j-Nb RHS j-Nb + F i pij + F j pjj + F k pjk; end if k > Nb RHS k-Nb RHS k-Nb + F i pik + F j pjk + F k pkk; end end plot the sparse matrix more exactly its sign then it is easier to see the pattern figure 6 ; imagesc sign abs A ; colormap 1-gray ; axis ij; axis equal; axis off; save as eps and svg needs the plot2svg function saveas gcf 'sparse_dirichlet eps' 'psc2' plot2svg 'Finite_element_sparse_matrix svg' ; calculate U then add zeros for the points on the boundary U_calc A\RHS; U_calc 0 1 Nb U_calc; exact solution u inline '1-x 2-y 2' 'x' 'y' ; U_exact u P 1 P 2 ; plot the computed solution figure 2 ; clf; plot_solution T P U_calc lw green save as eps and svg needs the plot2svg function saveas gcf 'computed_dirichlet eps' 'psc2' plot2svg 'Finite_element_solution svg' ; plot the exact solution figure 3 ; clf; plot_solution T P U_exact title 'Exact Dirichlet solution' saveas gcf 'exact_dirichlet eps' 'psc2' disp sprintf 'Error between exact solution and computed solution is 0 9g' max abs U_exact-U_calc given the l-th triangle calculate the integrals pij int_K lambda_i lambda_j and gij int_K grad lambda_i dot grad lambda_j function pii pij pik pjj pjk pkk gii gij gik gjj gjk gkk calc_elems i j k P1 P2 P3 extract the x and y coordinates of the points x1 P1 1 ; y1 P1 2 ; x2 P2 1 ; y2 P2 2 ; x3 P3 1 ; y3 P3 2 ; triangle area AT abs x2-x1 y3-y1 - x3-x1 y2-y1 /2; the integrals of lambda_i lambda_j pii AT/6; pij AT/12; pik AT/12; pjj AT/6; pjk AT/12; pkk AT/6; the integrals of grad lambda_i dot grad lambda_j gii x2-x3 x2-x3 + y2-y3 y2-y3 / 4 AT ; gjj x1-x3 x1-x3 + y1-y3 y1-y3 / 4 AT ; gkk x1-x2 x1-x2 + y1-y2 y1-y2 / 4 AT ; gij - x3-x1 x3-x2 + y3-y1 y3-y2 / 4 AT ; gik - x2-x1 x2-x3 + y2-y1 y2-y3 / 4 AT ; gjk - x1-x2 x1-x3 + y1-y2 y1-y3 / 4 AT ; function plot_solution T P U lw color fs Np k size P ; Nt k size T ; create a path in 3D that will trace the outline of the solution X zeros 5 Nt 1 ; Y zeros 5 Nt 1 ; Z zeros 5 Nt 1 ; for l 1 Nt i T l 1 ; j T l 2 ; k T l 3 ; m 5 l-4; X m P i 1 ; Y m P i 2 ; Z m U i ; X m+1 P j 1 ; Y m+1 P j 2 ; Z m+1 U j ; X m+2 P k 1 ; Y m+2 P k 2 ; Z m+2 U k ; X m+3 P i 1 ; Y m+3 P i 2 ; Z m+3 U i ; X m+4 NaN; Y m+4 NaN; Z m+4 NaN; end plot the solution plot3 X Y Z 'linewidth' lw 'color' color set gca 'fontsize' 15 function P get_points dummy_arg P 1 0 0 97799999999999998 0 20799999999999999 0 91400000000000003 0 40699999999999997 0 80900000000000005 0 58799999999999997 0 66900000000000004 0 74299999999999999 0 5 0 86599999999999999 0 309 0 95099999999999996 0 105 0 995 -0 105 0 995 -0 309 0 95099999999999996 -0 5 0 86599999999999999 -0 66900000000000004 0 74299999999999999 -0 80900000000000005 0 58799999999999997 -0 91400000000000003 0 40699999999999997 -0 97799999999999998 0 20799999999999999 -1 5 6700000000000004e-16 -0 97799999999999998 -0 20799999999999999 -0 91400000000000003 -0 40699999999999997 -0 80900000000000005 -0 58799999999999997 -0 66900000000000004 -0 74299999999999999 -0 5 -0 86599999999999999 -0 309 -0 95099999999999996 -0 105 -0 995 0 105 -0 995 0 309 -0 95099999999999996 0 5 -0 86599999999999999 0 66900000000000004 -0 74299999999999999 0 80900000000000005 -0 58799999999999997 0 91400000000000003 -0 40699999999999997 0 97799999999999998 -0 20799999999999999 2 6291902682773483e-17 -0 00065384615384583118 -0 29530436727611131 -0 40714990300538867 0 20468717553848803 -0 45950882973942581 0 5000339709144086 -0 052282439231331621 0 2044416534667938 0 45895712720185405 -0 37389783786652392 0 33573030516976354 -0 56520453690073769 -0 059175479864500904 0 55140920581395014 0 3176130751959379 -0 13138695201728712 0 62243041389833154 -0 1444609125190702 -0 68304604895205301 0 51896266786078527 -0 46675660322909635 -0 67677962425194449 0 22091656257348974 -0 61723328415727241 -0 35579831401388728 0 15396059620050356 -0 7270900369296075 0 73917109220890242 0 07757578859901848 0 70634528325694623 -0 23042511622333953 0 15390386456341718 0 72682700843038894 0 43656290869075642 0 60123684202225869 -0 43638925515227728 0 60099824488402309 -0 44644144213994863 -0 6148097863548887 -0 78986314125549695 -0 082937447632792677 0 48306021940714738 -0 66512339089274719 -0 71711649590963977 0 4137416136492385 0 7196608156574994 0 41521760024330068 -0 17239398182683685 0 81255391574260738 0 83205471084583948 -0 087400017493309945 0 3403176445451454 0 76427847186026787 -0 55164668399041139 0 34772453168234385 -0 58119856273351911 0 5140347159468972 -4 3977950621443274e-18 -0 83738246960490337 0 72531999275454517 -0 41850054828302335 -0 34213203269662834 -0 76835550876536474 -0 62483523734501689 -0 56238344018259578 0 56491731978950166 0 48814678592215621 0 38917253432269977 0 4180653756110998 0 23717545979432902 0 16912185462803336 -0 02596228168575496 0 28616935365733154 0 3938288333314971 0 24434831196145793 0 5426706476687696 0 13031124265795693 0 38005096874664096 0 07200907885348734 0 24522905811277632 -0 17819696266523521 0 22329825233727255 -0 0049527105560499022 -0 018125285606249857 -0 28382039345137922 -0 25412870038086932 -0 1265506189486573 -0 22001244253386174 0 12401796918826143 0 70594582936521699 0 24504049798687536 -0 28923932941907254 0 68786790036032985 -0 27673116940573905 0 49861546727758305 -0 19580411520445754 0 33991543187901008 -0 06480669781759768 0 45864775798987351 0 039811209497243899 0 59660089503236968 -0 012452577623367307 0 7574818502043279 0 36004239628137363 -0 5404259782287143 0 3198317338457855 -0 70312416491673124 0 2617529370159869 0 31866372248564073 0 11364902800796572 0 34001459209079515 -0 80612573707151369 -0 26251531242500947 -0 84358781498716018 0 088619865046719065 -0 44904385214632619 0 45802230787668152 0 35347968445219013 -0 076483406156372766 0 43529169673866364 -0 27201412098929839 0 33966312039894614 -0 38995864108359496 -0 44160050993872146 0 12473708976956265 0 22042251667383625 -0 60109068660701248 0 01773990267719644 -0 55195139956951267 0 57741967278428696 0 63172121266253345 0 83616296331029605 0 18384745602550973 -0 42194773652236522 0 74460017792571787 -0 15191200145175998 -0 50793571011947636 -0 28468886832575607 -0 60129916940148143 -0 15211412155416426 -0 35582508357904519 0 091328532274082594 0 85396134404515744 0 25268219760086502 -0 82436443072803334 -0 56522176794370282 0 65147861937265739 -0 17960003622540516 -0 84596380431778762 0 64262958611925292 -0 57845575520448644 -0 82977670051153962 0 27012165242582176 -0 75613161550480157 -0 43637469407737101 0 84682787874323973 -0 27560544844003693 0 66460522715024917 -0 070226896389436425 0 27135187110185849 0 61028308052190072 -0 5291687129171514 -0 72847571124389077 -0 46589500975993081 -0 44182174455360729 -0 44717423342713147 -0 23275286852251509 -0 91658877313709719 -0 096341120235654235 0 071238170548373975 0 15038100559175938 -0 07537518167909911 0 14811510831848701 0 14411045609387388 -0 30720135657166003 0 046251548820153726 -0 41167351754462878 0 080240926729569728 -0 1479533766093776 -0 090696855217268765 -0 13701177299039563 -0 14601986575596898 -0 0018723983808613188 -0 65587042533541751 -0 19614678101972682 -0 68222829844326438 0 056477001652015607 0 5820334614227094 -0 32431245274666742 0 56939947784766054 -0 17974699975143729 -0 33054333836607347 0 20546019299713006 -0 54252415463210835 0 21357850028719924 -0 397822414790862 -0 038177513011169929 -0 0015128273166917471 -0 6940216317375888 0 12814058587287039 -0 86580177941165415 -0 30513569295646459 -0 26238649633237193; function T get_triangles dummy_arg T 51 87 123 3 54 76 24 60 23 20 63 19 21 112 20 36 89 58 105 22 23 37 123 114 108 18 19 20 112 63 44 94 95 25 103 131 30 109 29 115 88 16 62 50 112 88 15 16 14 15 107 53 13 14 27 106 52 47 81 111 97 1 2 17 115 16 13 53 59 118 73 119 104 98 11 75 93 129 52 84 26 69 110 45 27 52 26 34 90 91 29 109 61 28 106 27 56 30 1 24 131 60 33 94 83 104 11 12 79 80 78 38 68 69 95 40 130 1 45 56 14 107 53 48 6 57 7 8 102 82 102 9 68 65 85 9 10 55 6 96 5 96 6 48 87 51 115 65 48 111 101 74 132 4 54 3 76 97 3 117 116 67 38 65 68 100 50 62 124 93 128 90 70 72 126 91 125 29 61 28 82 9 55 47 82 81 6 7 57 48 65 64 77 49 78 49 89 78 21 22 62 50 100 113 17 18 87 37 93 124 83 52 41 52 83 84 58 53 42 53 58 59 5 96 4 76 54 38 11 98 10 39 81 82 110 56 45 109 56 46 7 47 57 111 57 47 128 58 42 36 78 89 49 104 59 59 104 13 130 44 95 105 60 40 125 61 46 106 61 41 105 62 22 100 62 40 113 63 50 108 63 43 38 54 64 64 54 4 85 86 66 38 64 65 92 33 83 68 66 70 80 39 78 74 121 122 85 65 35 68 70 69 110 69 34 76 69 45 72 71 90 34 69 70 72 70 66 118 120 73 116 72 66 120 72 31 101 121 74 100 32 113 117 31 116 114 43 113 75 117 79 93 75 127 69 76 38 97 76 45 39 55 77 77 55 10 79 78 36 39 77 78 79 36 127 79 67 80 86 80 67 35 111 81 39 80 81 81 80 86 102 82 47 39 82 55 92 41 91 83 94 84 44 103 84 84 103 26 86 85 35 66 68 85 86 35 81 86 67 116 108 87 18 123 87 43 124 51 37 107 88 42 49 59 89 58 89 59 91 90 71 70 90 34 71 118 92 91 41 125 41 92 83 71 92 91 36 58 128 129 93 37 94 33 95 44 84 94 119 33 118 119 99 95 48 64 96 4 96 64 1 97 45 3 97 2 49 77 98 10 98 77 40 95 99 119 101 99 40 99 100 32 100 99 101 73 121 32 99 101 7 102 47 9 102 8 44 60 131 26 103 25 98 104 49 13 104 12 60 105 23 62 105 40 61 106 28 52 106 41 88 107 15 53 107 42 63 108 19 87 108 43 56 109 30 61 109 46 126 110 34 56 110 46 65 111 35 57 111 48 62 112 21 63 112 50 114 32 132 63 113 43 32 114 113 129 114 74 87 115 17 88 115 51 72 116 31 86 116 66 79 117 67 117 75 122 92 118 33 120 118 71 101 119 73 33 119 95 72 120 71 121 120 31 120 121 73 122 121 31 117 122 31 122 75 129 114 123 43 51 123 37 42 88 124 51 124 88 61 125 41 46 110 126 91 126 34 46 126 125 79 127 75 128 127 36 124 128 42 128 93 127 114 129 37 129 74 122 40 60 130 44 130 60 25 131 24 44 131 103 101 132 32 74 114 132; </source> Files by User Oleg Alexandrov from en wikipedia Triangulated irregular network Images with Matlab source code Finite element analysis | ||||
