Hi,
there is also a second system if you have constraints such that you can attach the new rows
to DOFs of the original matrix, namely A % ConstraintMatrix. An example of using
this to constraint values on a boundary to given function would be:
p => Solver % Variable % Perm
CM => AllocateMatrix()
CM % Format = MATRIX_LIST
ALLOCATE(temp(A % NumberOfRows), Indexes(Solver % MaxElementDOFs)); temp=0
DO t=1,Mesh % NumberOfBoundaryElements
Element=>ActiveBoundaryElement(t)
BC=>GetBC()
val(1:n) = GetReal(BC,'some setting',Found)
IF(.NOT.Found) CYCLE
n = GetElementNOFNodes()
nd = GetElementDOFs(Indexes)
CALL GetMassMatrix(Element,n,nd,val,MASS,FORCE) ! the usual elementwise
! integration, your own routine
DO i=1,nd
DO j=1,nd
CALL AddtoMatrixElement(CM,p(Indexes(i)),p(Indexes(j),MASS(i,j))
END DO
temp(p(Indexes(i)))=temp(p(Indexes(i))) + FORCE(i)
END DO
END DO
CALL List_toCRSMatrix(CM)
ALLOCATE(CM % Rhs(CM % NumberOfRows)); CM % Rhs=Temp(1:CM % NumberOFRows)
ALLOCATE(CM % InvPerm(CM % NumberOfRows))
DO i=1,CM % NumberofRows
CM % InvPerm(i) = i
END DO
A % ConstraintMatrix => CM
Note that you don't need to set the transpose, this is done automatically. Note also, that the CM matrix would most probably have a large number of zero rows, and one could (should?) compact the matrix using the permutation "InvPerm", not done here though.
The value of this code vs the addmatrix stuff is that the global matrix rows are only
shared by partitions as needed, and could be more scalable.
The code could be simplified if one could use DefaultUpdateEqations and friends
to update the CM matrix, unfortunately this is not really possible at the moment.
Something for the future perhaps.
The code is also completely unstested....
Regards,Juha
Integral COnstraint for single Solver
-
- Posts: 196
- Joined: 29 Sep 2011, 12:25
- Antispam: Yes
Re: Integral COnstraint for single Solver
HI there,
i have started working on the integral constraint again and i was hoping that someone could dive into the INvPerm for me a little bit. As i understand, in the pseudo code above there are ase many rows as in the original matrix. I would like to create a constraint matrix with just one row now. So i only add to the (1,perm(indexes(j)) entry and what do i do with the invperm then?
best regards
Franz
i have started working on the integral constraint again and i was hoping that someone could dive into the INvPerm for me a little bit. As i understand, in the pseudo code above there are ase many rows as in the original matrix. I would like to create a constraint matrix with just one row now. So i only add to the (1,perm(indexes(j)) entry and what do i do with the invperm then?
best regards
Franz
-
- Posts: 196
- Joined: 29 Sep 2011, 12:25
- Antispam: Yes
Re: Integral COnstraint for single Solver
Okay also without using the invperm i managed to make this work sor a serial run. In parallel it doesnt do the same thing though. It runs through but somehow the integral condition has a wrong effect. so i add the entries for all partitions in row number one and i add to rhs entry number one. i am realizing an integral condition so the processors should add up the entries in the right manner. is that also true for the rhs?
best regards
Franz
best regards
Franz