S.O.S. Mathematics CyberBoard

Your Resource for mathematics help on the web!
It is currently Mon, 28 Jul 2014 17:25:24 UTC

All times are UTC [ DST ]




Post new topic Reply to topic  [ 8 posts ] 
Author Message
PostPosted: Tue, 12 Jul 2011 14:36:32 UTC 
Offline
Math Cadet

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
 Profile  
 
PostPosted: Tue, 2 Aug 2011 18:08:50 UTC 
Offline
Member of the 'S.O.S. Math' Hall of Fame
User avatar

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.

2) Post your attempt
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
 Profile  
 
PostPosted: Fri, 12 Aug 2011 17:19:38 UTC 
Offline
Math Cadet

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:
http://www.imageupload.org/?d=46B3038B1
With best regards
freequo


Top
 Profile  
 
PostPosted: Fri, 12 Aug 2011 19:27:53 UTC 
Offline
Site Admin
User avatar

Joined: Sat, 26 Apr 2003 22:14:40 UTC
Posts: 2116
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
 Profile  
 
PostPosted: Sat, 13 Aug 2011 10:03:29 UTC 
Offline
Math Cadet

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
 Profile  
 
PostPosted: Sat, 13 Aug 2011 15:25:39 UTC 
Online
Moderator
User avatar

Joined: Wed, 30 Mar 2005 04:25:14 UTC
Posts: 13971
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
 Profile  
 
PostPosted: Mon, 15 Aug 2011 10:10:18 UTC 
Offline
Math Cadet

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
 Profile  
 
PostPosted: Thu, 17 May 2012 07:55:17 UTC 
Offline
Math Cadet

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
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 8 posts ] 

All times are UTC [ DST ]


Who is online

Users browsing this forum: No registered users


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

Search for:
Jump to:  
Contact Us | S.O.S. Mathematics Homepage
Privacy Statement | Search the "old" CyberBoard

users online during the last hour
Powered by phpBB © 2001, 2005-2011 phpBB Group.
Copyright © 1999-2013 MathMedics, LLC. All rights reserved.
Math Medics, LLC. - P.O. Box 12395 - El Paso TX 79913 - USA