S.O.S. Mathematics CyberBoard

Your Resource for mathematics help on the web!
 It is currently Mon, 1 Sep 2014 12:12:48 UTC

 All times are UTC [ DST ]

 Page 1 of 1 [ 8 posts ]
 Print view Previous topic | Next topic
Author Message
 Post subject: Stress return mapping algorithmPosted: Tue, 12 Jul 2011 14:36:32 UTC

Joined: Tue, 12 Jul 2011 14:31:58 UTC
Posts: 6
Dear all,
Firstly, I'm a new member, so if this problem is not posted in the right part of forum I hope that it will only be moved to right part and not deleted.

I'm having some problem with programming loop in mathematica. I'm working on an plasticity algorithm which has to be updated for any new step.

My opinion is that there should be at least three loops. The problem is that I don't know how to connect it.
I've manage to do first part of the loop, until F is less than 0. Loop also calculates next step when F becomes larger than 0. In that moment correction loop is needed.

Any help is welcome
Best regards
freequo

Top

 Post subject: Re: Stress return mapping algorithmPosted: Tue, 2 Aug 2011 18:08:50 UTC
 Member of the 'S.O.S. Math' Hall of Fame

Joined: Thu, 7 Apr 2005 16:32:15 UTC
Posts: 1276
Location: United States EST
freequo wrote:
I'm having some problem with programming loop in mathematica. I'm working on an plasticity algorithm which has to be updated for any new step.

My opinion is that there should be at least three loops. The problem is that I don't know how to connect it.
I've manage to do first part of the loop, until F is less than 0. Loop also calculates next step when F becomes larger than 0. In that moment correction loop is needed.

freequo, Welcome to sosmath!

You'll probably have to be a bit more specific. I've tried a search for "plasticity algorithm" and got quite a diverse array of options. Here are some suggestions if you'd like some help.

1) Post the problem you're trying to program.

You can post code by using the 'code' tags:

Code:
#include<iostream>
using namespace std;

int main
{
cout << "Hello World! \n";
return 0;
}

Hit quote to see what I typed in.

Hope this helps out

Best,

-Qleak

Top

 Post subject: Re: Stress return mapping algorithmPosted: Fri, 12 Aug 2011 17:19:38 UTC

Joined: Tue, 12 Jul 2011 14:31:58 UTC
Posts: 6
Yeah you are wright.
I thought that at least I put an image of an flow chart.
Here is a code written in Mathematica:
Code:
state[\[Sigma]10_, \[Sigma]20_, e0_, n0_, c0_, f0_, step_] :=
Module[{matD, F, c, \[Phi], \[Lambda], I1, J2, \[Sigma]1, \[Sigma]2,
Em, \[Nu], matDep, \[Sigma], edf, avek},
c = c0; \[Phi] = f0; \[Sigma]1 = \[Sigma]10; \[Sigma]2 = \[Sigma]20;
Em = e0; \[Nu] = n0;

edf = Table[{If[i == 1, 0, 0.00001], 0}, {i, 1, step}];

\[Sigma] = {{\[Sigma]1}, {\[Sigma]2}};

matD = Em/((1 + \[Nu]) (1 - 2 \[Nu])) ({
{1 - \[Nu], \[Nu]},
{\[Nu], 1 - \[Nu]}
});

Do[\[Sigma] = \[Sigma] + edf[[i]].matD;

F = -((6 c Cos[\[Phi]])/(Sqrt[3] (Sin[\[Phi]] + 3))) + Sqrt[
1/6 ((\[Sigma][[1, 1]] - \[Sigma][[2, 1]])^2 + \[Sigma][[1,
1]]^2 + \[Sigma][[2, 1]]^2)] + ((\[Sigma][[1,
1]] + \[Sigma][[2, 1]]) (2 Sin[\[Phi]]))/(
Sqrt[3] (Sin[\[Phi]] + 3));

term = If[F <= 0, "Elastic", "Plastic"];

avek = If[\[Sigma][[1, 1]] == \[Sigma][[2, 1]] == 0, 0, ({
{(2 \[Sigma][[1, 1]] +
2 (\[Sigma][[1, 1]] - \[Sigma][[2, 1]]))/(
2 Sqrt[6]
Sqrt[\[Sigma][[1,
1]]^2 + (\[Sigma][[1, 1]] - \[Sigma][[2,
1]])^2 + \[Sigma][[2, 1]]^2]) + (2 Sin[\[Phi]])/(
Sqrt[3] (3 + Sin[\[Phi]]))},
{(-2 (\[Sigma][[1, 1]] - \[Sigma][[2, 1]]) +
2 \[Sigma][[2, 1]])/(
2 Sqrt[6]
Sqrt[\[Sigma][[1,
1]]^2 + (\[Sigma][[1, 1]] - \[Sigma][[2,
1]])^2 + \[Sigma][[2, 1]]^2]) + (2 Sin[\[Phi]])/(
Sqrt[3] (3 + Sin[\[Phi]]))}
})];

\[Lambda] = If[F > 0, F/Flatten[Transpose[avek].matD.avek][[1]], 0];

matDep = If[F > 0, matD (( {
{1, 0},
{0, 1}
} ) - avek.Transpose[avek].matD/
Flatten[Transpose[avek].matD.avek][[1]]), matD];

Print[{\[Sigma] // MatrixForm, F, term, \[Lambda],
matDep // MatrixForm}], {i, 1, step}];

]

And here is an image of flow chart:
With best regards
freequo

Top

 Post subject: Re: Stress return mapping algorithmPosted: Fri, 12 Aug 2011 19:27:53 UTC

Joined: Sat, 26 Apr 2003 22:14:40 UTC
Posts: 2119
Location: El Paso TX (USA)
Mathematica is a function-oriented language. Whenever possible you should avoid loops.

_________________
The greater danger for most of us lies not in setting our aim too high and falling short; but in setting our aim too low, and achieving our mark. - Michelangelo Buonarroti

Top

 Post subject: Re: Stress return mapping algorithmPosted: Sat, 13 Aug 2011 10:03:29 UTC

Joined: Tue, 12 Jul 2011 14:31:58 UTC
Posts: 6
Didn't know that. Thanks for the advice.
I don't know to use any other problem, and I need this loop for programming finite element.
Thanks

Top

 Post subject: Re: Stress return mapping algorithmPosted: Sat, 13 Aug 2011 15:25:39 UTC
 Moderator

Joined: Wed, 30 Mar 2005 04:25:14 UTC
Posts: 14005
Location: Austin, TX
freequo wrote:
Didn't know that. Thanks for the advice.
I don't know to use any other problem, and I need this loop for programming finite element.
Thanks

That really depends, a lot of times you can avoid doing loops in mathematica by clever use of their functions. I remember when I was first starting out on it I did lots of loop-y things too, but every time I managed to find a better, cleaner way to do it w/o loops.

_________________
(\ /)
(O.o)
(> <)
This is Bunny. Copy Bunny into your signature to help him on his way to world domination

Top

 Post subject: Re: Stress return mapping algorithmPosted: Mon, 15 Aug 2011 10:10:18 UTC

Joined: Tue, 12 Jul 2011 14:31:58 UTC
Posts: 6
I've made some program for some other geotechnical problem using functions, and only functions, but I think that this problem could not be solved without looping.
Anyway thank you for the advice.
freequo

Top

 Post subject: Re: Stress return mapping algorithmPosted: Thu, 17 May 2012 07:55:17 UTC

Joined: Tue, 12 Jul 2011 14:31:58 UTC
Posts: 6
Here it is:

Code:
return[tol_, sini_, \[Sigma]y_, Em_, nu_, status_] :=
Module[{matD, \[Nu] = nu, \[Lambda], s = sini, matDep, cep},

F[x_, y_, z_] := \[Sqrt]((x)^2 + (y)^2 - x y + 3 (z)^2) - \[Sigma]y;

avek[x_, y_, z_] := ({
{(2 x - y)/(2 Sqrt[(x)^2 - x y + (y)^2 + 3 (z)^2])},
{(-x + 2 y)/(2 Sqrt[(x)^2 - x y + (y)^2 + 3 (z)^2])},
{(3 z)/Sqrt[(x)^2 - x y + (y)^2 + 3 (z)^2]}
});

matD = Em/(1 - \[Nu]^2) ({
{1, \[Nu], 0},
{\[Nu], 1, 0},
{0, 0, 1 - \[Nu]}
});
matDep = matD (( {
{1, 0, 0},
{0, 1, 0},
{0, 0, 1}
} ) -
Flatten[Transpose[avek[s[[1]], s[[2]], s[[3]]]].avek[s[[1]],
s[[2]], s[[3]]]][[1]] matD/
Flatten[Transpose[avek[s[[1]], s[[2]], s[[3]]]].matD.avek[
s[[1]], s[[2]], s[[3]]]][[1]]);

cep = If[status === "Elastic", matD, matDep];

While[F[s[[1]], s[[2]], s[[3]]] > tol,
avek[s[[1]], s[[2]], s[[3]]];
\[Lambda] =
Flatten[F[s[[1]], s[[2]], s[[3]]]/
Flatten[Transpose[avek[s[[1]], s[[2]], s[[3]]]].cep].avek[
s[[1]], s[[2]], s[[3]]]][[1]];

s = Flatten[s - \[Lambda] cep.avek[s[[1]], s[[2]], s[[3]]]];

Print[F[s[[1]], s[[2]], s[[3]]]];

];

{s, av = avek[s[[1]], s[[2]], s[[3]]]}

]

Top

 Display posts from previous: All posts1 day7 days2 weeks1 month3 months6 months1 year Sort by AuthorPost timeSubject AscendingDescending
 Page 1 of 1 [ 8 posts ]

 All times are UTC [ DST ]

Who is online

Users browsing this forum: No registered users

 You cannot post new topics in this forumYou cannot reply to topics in this forumYou cannot edit your posts in this forumYou cannot delete your posts in this forum

Search for:
 Jump to:  Select a forum ------------------ High School and College Mathematics    Algebra    Geometry and Trigonometry    Calculus    Matrix Algebra    Differential Equations    Probability and Statistics    Proposed Problems Applications    Physics, Chemistry, Engineering, etc.    Computer Science    Math for Business and Economics Advanced Mathematics    Foundations    Algebra and Number Theory    Analysis and Topology    Applied Mathematics    Other Topics in Advanced Mathematics Other Topics    Administrator Announcements    Comments and Suggestions for S.O.S. Math    Posting Math Formulas with LaTeX    Miscellaneous