changing bc during solution

[Ali Azary]

hi
I need to change boundary conditions during solution. as stresses pass a critical value at nodes I want to set the boundary and displacement there and then solve for the new boundary the second part. I wrote a user macro and edited reference arrays starting at  mr(np(31)+nneq) and hr(np(27)+ndf*numnp). plotting the boundary shows that the boundary conditions are properly changed. However it seems that to make the program to enforce the new bc and solve for it needs more than that. I examined the codes in routine pmacr2.f, and the ones called therein such as pmesh.f and pbouin.f . The codes are extremely complicated and I was not able to get anywhere with them. I would appreciate it if you could help me with the code needed to enforce  new bc or any changes made in main arrays in general.

thank you for your time and consideration.

[Prof. S. Govindjee]

After you change the arrays you have to tell FEAP to recompute the equation numbers etc., followed by re-equilibrating the state of the system.
In PLOT there is a method for doing this interactively (plot> pbou) see the manual for how to use the feature.  Also look at plot/pplotf.f and search for pbou to see how to tell FEAP to reset the equation numbers.

[Prof. R.L. Taylor]

There is a simple way to compute the new numbers.

In your routine, add the common block region using
              include 'region.h'

              nreg = -1 ! forces all regions in analysis to be renumbered

              call pnewpr() ! renumbers all equations and resets profile

You should then have the new numbering for the next solution.

[Ali Azary]

thank you very much for your help. I tried what you said that I realized it is the first thing that is done in pmacr3.f. It does not work properly. I attached my macro and input files for your consideration. I had wrote the new bc in a file and used the "mesh file" command as you can see. since the extra read write slows down the code execution, I wanted to implement it in the code directly. I am so grateful for your kind help. I will try to dig more into the problem with the kind instructions you and dr Govindjee provided. however I would appreciate it if you could take a look at my files and guide me in the right direction.

thank you very much again for your time and kind consideration

[Ali Azary]

I just had a look at pnewpr.f and in the beginning there was a call to profil(mr(np(20+k1), ... for k1=1,4 . In palloc.f for pointer 21 through 24 said profile pointer to partitions 1 to 4. I thought this might point to the dof partitions used for staggered solution scheme. so since I just change the bc for the second partition I thought I just need to reset the data for the second partition. So , i copied the subroutine pnewpr and appended that to my macro and changed the name. I just changed the loop "do k1=1,4" to "k1=2,2". However, it did not work. I would appreciate it if you could tell me what is wrong. It probably has something to do with proportional loading for displacement bc for partition 1. when calling pnewprf for all partitions it seems to neglect increasing bc with time and enforces it all at once.

Thank you very much

[Prof. R.L. Taylor]

What version of feap are you using?  The date on the umacr is 2004.

It is not clear what you want to do when you reset the equations.  The current crack routine has all the computations commented and you just set all the dof=3 equations to "fixed" and the value to 1.  So that does not seem right.  You should do a very small problem in which you use crack.  Before crack use the command 'show id' to output all the values in the boundary condition array; repeat after the 'crack' command. Do the numbers become what you want.  Did the boundary condition change and were the new numbers correct.  Once that works you should be able to chack the rest.

If this is a very old version of feap it may be that there were bugs in changing boundary conditions during analysis; in the current release it may be necessary to use an 'opti' command after each boundary condition change to fore the pnewpr call.

Of minor importance now, your input file should always have a 'stop' command at the location where the analysis is complete.

[Ali Azary]

Dear Dr Taylor

Thank you very much for your reply. Our feap version is 8.3. I hope it is not a bug related to this version. I am sorry if my explanations was not clear enough and the files were quite messy. I tried "optimize" command after "crack" command that changes BC as you instructed and it seems to do the trick. However there is one minor issue come up and probably is a bug in the code: in contour plots it seems that a node is being blown out, but seemingly it has no effect on the results. I attached the images for your consideration. Regarding my problem, I am solving cracking problem using phase-field method, which is basically gradient damage in fact. the first two DOFs are displacements and the third DOF is the damage variable. I first solve the elasticity problem and determine the stresses, then I want to check if the maximum tensile stress is passed the critical value. if so the point is cracked so I set the damage value to one, as BC for the third DOF. Then I solve the damage equation to determine the regularized crack. I hope I have been clear enough, and please excuse me if I have been verbose. In case you would like to take a look I attached the new code and input file, too. I cleaned up the code and added explanatory comments. Thank you again for your kind help.

[Prof. R.L. Taylor]

Can you display allocation line for ID form pcontr.f  - how big is the allocation?

Can you look in the output file and see what the code says for partitions.

Finally, can you tell what node is causing the spurious display?  Or is it a deformed plot using some value for the displacement (maybe the damage variable?)  The command PLOT DOFS 1 2 0 would mask off the damage variable.

[Ali Azary]

the allocation is as below:
l2   = max(ndf*numnp,1)
setvar = palloc( 31,'ID   ',l2*2        ,  1)

[Ali Azary]

It happened in one case, other trials did not show the strange behavior. However I realized the results are way different from when I used to write the BC in a file in the user macro and use mesh file command afterward. The output file seems normal to me. I am not sure what you are looking for, so I attached an output n case you want to take a look. I remember I checked the displacements but no displacement seemed so off. Regarding the plots they are undeformed. I used the input file "iplate" I attached in my last message. I appreciate your help and I will understand if you do not want to put more time. You have been very helpful so far and I am so grateful for that.

[Ali Azary]

I keep getting server errors. I broke the message into parts. However, i could not attach the output file.

[Prof. R.L. Taylor]

In checking the output file just make sure that the partition command is correctly displayed in the output.  The manner in specifying partitions changed during v8.3 and also the storage for the ID array changed.  There may be more than one allocation location for the ID array: the initial one when all other arrays are allocated and a second one after partitions are allocated.  In particular we now store a copy of the ID array for each partition.  I just want to know how your program is storing the array.  One way to discover is to enter the program interactively after all the data is given and before you solve any problems.  Use the command SHOW ID to display the array.  The program should tell you how many parts exist for your problem (normally this is only two, but with partitions it is now more than two).  My version of 8.3 has an allocation in pcontr.f at line 817 that allocates npart+1 copies of the ID array.  The storage in the array is a bit awkward as it was changed form an earlier form:
(1) the first part id(ndf,numnp,1) stores the equations for partition 1 (this is the storage when no staggers are used)
(2) the second part id(ndf,numnp,2) stores the boundary codes for all the dofs, a 0  or an active dof and a non-zero for a fixed or unused one
(3) the remaining parts id(ndf,numnp,3) etc. store the equation numbers for partiaion 2 etc.
You should be able to follow how this is set up in the lines after the allocation in pcontr.f

So if the equations are getting set correctly (you can break the batch solution after every CRACk and OPTI command) and use an interactive solution command SHOW ID to look - then at least you know you are setting the equations correctly.

The remainder is to set the values for the solution. It appears you are setting these in the array "F" np(27) with an offset; so SHOW F should tell you if they are correctly done.

Once this seems correct then check if the solution is as expected.  The process you are doing is fairly complex and it is possible that some other changes are still needed -- but the above is what is done when we use MESH commands to change BOUN and DISP values (or their E or C alternatives).  I have noted that the current version of feap does not force a renumber of equations if a BOUN type change is made.  An improvement will be implemented in future releases.

[Ali Azary]

thank you very much. I will try it out, and I think I might be able to figure out how it works with the kind instructions you provided.