I'm a composite materials person so come from a different perspective (although bone is essentially a composite of HCA and collagen!). Doing this by (ab)using the .rpy scripts created by Abaqus sounds like a nightmare to me and I can't think why your supervisor would specify this. However, you could write python programs to do what you are after though (or MatlLab, we use both for this at my company)
Anyway, AFAIK Abaqus always requires a section to be assigned to individual elements in order for these elements to each have their own properties assigned to them, even if you're using continuum elements. As a result, if you want to update/assign properties to individual elements than each element must have its own section assigned to it.
Creating a section for every element who's properties you want to change via CAE or from .rpy files would be extremely tedious. If you can get the element numbers of all the elements you wish to modify then I think this should be quite doable . As a start you could just pick the elements from CAE with some jedi mouse clicking skills or get creative with partitions.
Once you know the numbers of elements of interest you could then write a script to create a material (and an orientation, seeing as bone is a composite) for every element/section, assign them and then get the program to print and slice the relevant lines into your input deck.
All the values for the materials could start off a being identical. What you could then do is run the model and use the script you have already created to identify elementsthat have exceeded your threshold parameter. You would then get your program from above to change the material properties associated with those elements. The re-run it.
As an example; for creating materials, pseudo code might be something like:
Nel= number of elements you are modifying
Create array or list or something of the element numbers of interest
Map, index, or similar (whatever is your preference) element numbers to list of 1 to Nel
Then start a loop such as;
while x is less than or equal to Nel:
print '*Material, name=Material-',x
print '*Elastic, type=ENGINEERING CONSTANTS'
print ' Elastic constant xA, Elastic constant xB, Elastic constant xC, etc...'
increase x by one
Reset x once you're out the loop
Map output back to element numbering list
Slice back into Abaqus .inp file
Once you've got something along the above lines running then write a program to do this automatically by going back and forth from Abaqus for a given number of iterations.
If you want to improve on what you've outlined (IMO at least) then what you could do is get Abaqus to print out the stress components for each element, read these from the .odb (there are various post on the internet on how to do this), take their eigen values and vectors, update material properties in an anisotropic manner, as that is what bone would naturally do.
