SCOREC / fep

Finite Element Programming course materials
6 stars 4 forks source link

integration point fields #11

Closed cwsmith closed 3 years ago

cwsmith commented 5 years ago

@RuixiongHu writes:

  1. I tried to initialize a field on integration point, I did it. However, in postprocessing I figured out this "state variable" value I set it to be is only correctly shown in 'x' component. But 'magnitude' is somewhat little different, seems has been adapted according to coordinate value. I initialize it by following call, there's only one component and it's a scalar, and I printed out my assigned value ‘stateptr’, which is correct. So I was wondering why what I entered is not into 'Magnitude' but instead was shown by 'x' component.
  pField SV = pumi_field_create( mesh, "StateVariable", 1,1, pumi_shape_getIP(3,1) ); 
  pumi_node_setField( SV, e, 0, stateptr );
  1. Can I know how to iterate over integration point? I don't see any example related in the documentation. I was using a tetrahedral mesh with one integration point. (Initialized by above)

  2. Currently what I have is a simple 3D block with tetrahedral mesh. Is it possible for me to use pumi_mesh_createElem( ) to create some 2D surface element to it? I could delete them when writing output VTK. What I would like to do is I have a few point, I want to connect them and form a surface. So I'm creating element and save those 2d element into a vector so that I could loop over them later. PS: the 3rd question, if I create 2D element into this 3D model, it gives me error signal 6 (abort3) during calling mpi_mesh_freeze(), probably because this is a 3D mesh file. So instead, inside this current file, can I create a new pMesh object and save the new 2D mesh into it? Just like the construction example, the only difference is I have two pMesh objects, one is a 3D mesh file and another is a 2D mesh file.

RuixiongHu commented 5 years ago

Hi just wanna do a quick update, if I create a new pMesh object. Then when I try to iterate over one of them, it give me:

i < m->end[t] failed at /gpfs/u/barn/CCNI/CCNIsmth/systemSoftware/core/mds/mds.c + 590

Thanks a lot!

mortezah commented 5 years ago

As far as I can tell this statement pField SV = pumi_field_create( mesh, "StateVariable", 1,1, pumi_shape_getIP(3,1) ) creates a vector field with one component.

If you need to create a scalar field use pField SV = pumi_field_create( mesh, "StateVariable", 1,0, pumi_shape_getIP(3,1) )

mortezah commented 5 years ago

You could use this to get the field value https://github.com/SCOREC/core/blob/4446cdaeafa2f86b11444e8d256ac4d91bb25d18/pumi/pumi_field.cc#L171

So for getting the integration point you would iterate over tets and use the above function to get the value at each node (which in this case would be the integration point)

mortezah commented 5 years ago

In regards to your 3rd question: I really don't understand what is it you are trying to achieve. In theory, you should be able to add new entities of any type to the mesh but you have to note that that may make your mesh invalid (meaning that the verification steps may not pass).

RuixiongHu commented 5 years ago

@mortezah

Hi thanks a lot for the explanation! The field create is successful now.

And for the loop over integration point, I see that function could get me the value of integration point. But can I also ask for is there functions that I could get the location (coordinate value) of the integration point?

For the 3rd question. Yes I tried, I could add new entities of any type and the validation steps are giving error. I'm working in a 3D mesh, but wanna add some 2d surface elements into it, so that inside the code I could know their coordinate maybe later. For example I would like to know the (x,y,z) of this surface to update my field. They don't have to be in the output file. Is there any better way to do it?

Thanks, I'll really appreciate your help

mortezah commented 5 years ago

To get the location of the nodes in the parent entity (i.e., xi_0, xi_1, xi_2) we have the following in apf

virtual void apf::FieldShape::getNodeXi | (int type, int node, Vector3 & xi)

But that functionality does not seem to be exposed in pumi

Also, note that these are standard integration rules. For example, for a first order tet the integration point will be at the parametric location (0.25, 0.25. 0.25).

mortezah commented 5 years ago

Still, I am not clear what is needed in regards to your 3rd question. You said, "I'm working in a 3D mesh, but wanna add some 2d surface elements into it, so that inside the code I could know their coordinate maybe later." In order to create those entities, you need to know the coordinates anyway. Why don't you just keep track of those coordinates as opposed to making new entities as a way of tracking those coordinates?

RuixiongHu commented 5 years ago

Hi!

Thanks a lot for the help, I'll be thinking over to come up with a better way to implement this.

Can I ask a one last question?

Is there functions in pumi that could project the integration point field onto nodes? From other FEM solver I could get my variable of interest at integration point. And I would like to have those projected onto nodes for a smooth transition.

Thanks and look forward for the reply

mortezah commented 5 years ago

I am not aware of such a function in pumi. I should also say that I don't have much experience with pumi. I mostly use apf directly.