Change boundary condition during calculation

[greyybear]

Dear Professor Taylor,
dear FEAP Admin,

I am currently trying to couple FEAP with a flow solver and an optimization code at the same time.  For testing purposes I created the mesh of a half circle which is fully restraint at its center line.

While the flow solver calculates forces that will be applied at the exterior surface of the half circle, the optimization code will calculate new displacements of the exterior surface, which lead to a new interior mesh design. For the resulting new shape (e.g. a half ellipse) the flow solver will now calculate new forces and so on... As this is an iterative process, the boundary type of the exterior surface would have to be changed periodically from forces to displacements and vice versa.

I  tried to change the boundary type manually in a user macro but wasn't able to issue a remeshing/resetting of the problem. I would be really thankful if you could point me in the right direction.

Sincerely,

Mani Gebbers

[Prof. R.L. Taylor]

Changing boundary conditions involves more steps than you are probably performing.  We suggest you look at all the steps we perform when you call 'mesh' as a solution command.  All the steps can be found in the subprogram pmacr3.f (in the /program directory) under the comment that contains [mesh.

You need to have feap reset the equation numbers in paraticular -- some of the other steps may not be necessary.

You can also see that there is an undocumented option for a solution command
       MESH filename
I would expect you could create two files: one for force b.c. and the other for displacement ones.

You could then have your commands control whcih b.c. to use and feap will reset the equation numbers.  One way to watch what is happening is to issue the commands in an interactive mode
and look at the arrays for b.c./equations (ID array which is labeled as 31 if you use the SHOW DICT solution command).  This can also help you to debug your user routine by looking to see if you are changing things correctly.

Good luck.

[greyybear]

Dear Professor Taylor,

thank you very much for your quick reply! I had a look at the "MESH, filename" solution command because it seemed to be exactly what i was looking for.

I created two files:
1st file: contains force b.c and the solution command (ihalfcircle)
2nd file:contains displacement b.c. (idisp_bc)

When I use the solution command written below, the force boundary problem is solved. And as I added a INTEr command at the end it is still possible to enter additional solution commands. If now the "MESH, filename" command is entered, the second mesh file is read and the problem is set up properly as a displacement problem. Entering another tang,,1 solves the displacement boundary problem correctly.

Unfortunately, I wasn't able to use the "MESH, filename" command in a pure BATCh mode. So, when I remove the comment in the solution command below, the following tang,,1 command is ignored and no displacements for the inner grid nodes are calculated.  Am I making an apparent mistake?

Regards,
Mani Gebbers

Code: [Select]

BATCh
!MESH, idisp_bc
 LOOP,newton,10
   TANG,,1
 next,newton
 DISPlacement all
 plot,wipe
 plot,defo
 plot,mesh
 plot,cont,2
END
INTEr
STOP


[Prof. R.L. Taylor]

I put together an example that does your problem with first force b.c., then displ, then force again and then displ.

Note the use of the BOUNd SET - this forces b.c. to go back to zero for the next force solution.

[Prof. R.L. Taylor]

HEre is the file

[greyybear]

Dear Professor Taylor,

With the help of the input file you provided I was able to identify my mistake. I have successfully coupled FEAP with an optimization and flow solving code. For small to medium deformations everything works fine, unfortunately for large deformations the FEAP displacement problem does not always converge and thus leads to degenerated meshes.

That's why I switched to another way for updating the internal mesh in the end. Thank you very much for your help nonetheless, it is much appreciated.

Mani Gebbers

[Yijun Chen]

Dear Professors,

Thank you very much for posting the mesh.zip file, It works well using the feap command.
Then I modified the mesh.zip file, and try to use the parfeap command to solve the problem. I found the result seems not to work well.

I use version 8.3.

Could you point me in the right direction in your free time? Thank you very much.

Yijun Chen

[Prof. S. Govindjee]

Doing this in parfeap will be somewhat more tricky.  The added complication is that when you reset the boundary conditions you will need to recompute all the mesh partitioning, equation numbering, etc.  Right now this is done during the serial run of parfeap/feap.  There is also a big dependency on the if you keep the driven dofs in the system of equations or not.  If you keep them in, then perhaps you do not need to rerun the partitioning etc.  But I would have to think carefully about this.  If you are using something like mesh,filename also beware that you will need separate files for each partition of your problem.   It could get tricky to set up but seems like it is possible.

[Yijun Chen]

Dear Prof. Govindjee,

Thank you very much for your guidance.

I want to change the boundary condition from 1 to 0. The results seem not to reflect this change. I only use one partitioning for investigating this problem.
I will also think carefully about his problem.

Best regards,

Yijun Chen 

[Prof. S. Govindjee]

If you are doing a one processor job, I think your case may work 'out of the box' if you create your single partition input file using something like GRAPh,NODE,1  followed by OUTDomains,AIJ,1

[Yijun Chen]

Dear Prof. Govindjee,

Thank you very much for your reply.

I try this strategy again, but the result for the reset mesh also seems not correct. I post some files.
For the second mesh, I reset the x-direction of nodes 3, 6, and 9 into free. I guess the x-direction displacement of nodes 3, 6, and 9 should be changed. But based on the second mesh.png,
these values are still the same as that in the first mesh.png.
So for the parfeap command, I guess we cannot reset the boundary condition (degree of freedom) from 1 to 0 (constraint/prescribed displacement to free).

Best regards,

Yijun Chen

[Yijun Chen]

For the feap command framework, I think the results are no problem.
I also post these files.

[FEAP_Admin]

I see the problem.  Version 8.3 does not have the AIJ,1 option on OUTDomains.  I suggest upgrading to a current release of FEAP.  This is especially important if you want to use current PETSc releases.

One further suggestion that may work with 8.3 is to try OUTD,BLOC; in 8.3 the BLOCked equation option includes all dofs in the assembled matrix so it may be possible to still switch boundary conditions.  My  only concern is that I do not recall if we took this into account during the pre-allocation of memory.
Also, when you use BLOCked your matrix type will be BAIJ which does not work with all PETSc solvers.

Last comment, in the files you uploaded the serial problem is not the same as the parfeap problem.  I suggest using the same problem to make it easier to see if things are working correctly.

[Yijun Chen]

Dear Professor,

Thank you very much for your guidance.

I am going to upgrade to FEAP 8.6 after PhD graduation.

I try the OUTD,BLOC option, and the results also not correct.

Best regards,

Yijun Chen