1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/ 2 #include <petscdmplex.h> 3 #include <petscblaslapack.h> 4 5 PetscErrorCode DMPlexGetTransitiveClosure_Internal(DM, PetscInt, PetscInt, PetscBool, PetscInt *, PetscInt *[]); 6 7 struct _n_Petsc1DNodeFamily 8 { 9 PetscInt refct; 10 PetscDTNodeType nodeFamily; 11 PetscReal gaussJacobiExp; 12 PetscInt nComputed; 13 PetscReal **nodesets; 14 PetscBool endpoints; 15 }; 16 17 static PetscErrorCode Petsc1DNodeFamilyCreate(PetscDTNodeType family, PetscReal gaussJacobiExp, PetscBool endpoints, Petsc1DNodeFamily *nf) 18 { 19 Petsc1DNodeFamily f; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 ierr = PetscNew(&f);CHKERRQ(ierr); 24 switch (family) { 25 case PETSCDTNODES_GAUSSJACOBI: 26 case PETSCDTNODES_EQUISPACED: 27 f->nodeFamily = family; 28 break; 29 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family"); 30 } 31 f->endpoints = endpoints; 32 f->gaussJacobiExp = 0.; 33 if (family == PETSCDTNODES_GAUSSJACOBI) { 34 if (gaussJacobiExp <= -1.) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Gauss-Jacobi exponent must be > -1.\n"); 35 f->gaussJacobiExp = gaussJacobiExp; 36 } 37 f->refct = 1; 38 *nf = f; 39 PetscFunctionReturn(0); 40 } 41 42 static PetscErrorCode Petsc1DNodeFamilyReference(Petsc1DNodeFamily nf) 43 { 44 PetscFunctionBegin; 45 if (nf) nf->refct++; 46 PetscFunctionReturn(0); 47 } 48 49 static PetscErrorCode Petsc1DNodeFamilyDestroy(Petsc1DNodeFamily *nf) { 50 PetscInt i, nc; 51 PetscErrorCode ierr; 52 53 PetscFunctionBegin; 54 if (!(*nf)) PetscFunctionReturn(0); 55 if (--(*nf)->refct > 0) { 56 *nf = NULL; 57 PetscFunctionReturn(0); 58 } 59 nc = (*nf)->nComputed; 60 for (i = 0; i < nc; i++) { 61 ierr = PetscFree((*nf)->nodesets[i]);CHKERRQ(ierr); 62 } 63 ierr = PetscFree((*nf)->nodesets);CHKERRQ(ierr); 64 ierr = PetscFree(*nf);CHKERRQ(ierr); 65 *nf = NULL; 66 PetscFunctionReturn(0); 67 } 68 69 static PetscErrorCode Petsc1DNodeFamilyGetNodeSets(Petsc1DNodeFamily f, PetscInt degree, PetscReal ***nodesets) 70 { 71 PetscInt nc; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 nc = f->nComputed; 76 if (degree >= nc) { 77 PetscInt i, j; 78 PetscReal **new_nodesets; 79 PetscReal *w; 80 81 ierr = PetscMalloc1(degree + 1, &new_nodesets);CHKERRQ(ierr); 82 ierr = PetscArraycpy(new_nodesets, f->nodesets, nc);CHKERRQ(ierr); 83 ierr = PetscFree(f->nodesets);CHKERRQ(ierr); 84 f->nodesets = new_nodesets; 85 ierr = PetscMalloc1(degree + 1, &w);CHKERRQ(ierr); 86 for (i = nc; i < degree + 1; i++) { 87 ierr = PetscMalloc1(i + 1, &(f->nodesets[i]));CHKERRQ(ierr); 88 if (!i) { 89 f->nodesets[i][0] = 0.5; 90 } else { 91 switch (f->nodeFamily) { 92 case PETSCDTNODES_EQUISPACED: 93 if (f->endpoints) { 94 for (j = 0; j <= i; j++) f->nodesets[i][j] = (PetscReal) j / (PetscReal) i; 95 } else { 96 for (j = 0; j <= i; j++) f->nodesets[i][j] = ((PetscReal) j + 0.5) / ((PetscReal) i + 1.); 97 } 98 break; 99 case PETSCDTNODES_GAUSSJACOBI: 100 if (f->endpoints) { 101 ierr = PetscDTGaussLobattoJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);CHKERRQ(ierr); 102 } else { 103 ierr = PetscDTGaussJacobiQuadrature(i + 1, 0., 1., f->gaussJacobiExp, f->gaussJacobiExp, f->nodesets[i], w);CHKERRQ(ierr); 104 } 105 break; 106 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Unknown 1D node family"); 107 } 108 } 109 } 110 ierr = PetscFree(w);CHKERRQ(ierr); 111 f->nComputed = degree + 1; 112 } 113 *nodesets = f->nodesets; 114 PetscFunctionReturn(0); 115 } 116 117 static PetscErrorCode PetscNodeRecursive_Internal(PetscInt dim, PetscInt degree, PetscReal **nodesets, PetscInt tup[], PetscReal node[]) 118 { 119 PetscReal w; 120 PetscInt i, j; 121 PetscErrorCode ierr; 122 123 PetscFunctionBeginHot; 124 w = 0.; 125 if (dim == 1) { 126 node[0] = nodesets[degree][tup[0]]; 127 node[1] = nodesets[degree][tup[1]]; 128 } else { 129 for (i = 0; i < dim + 1; i++) node[i] = 0.; 130 for (i = 0; i < dim + 1; i++) { 131 PetscReal wi = nodesets[degree][degree-tup[i]]; 132 133 for (j = 0; j < dim+1; j++) tup[dim+1+j] = tup[j+(j>=i)]; 134 ierr = PetscNodeRecursive_Internal(dim-1,degree-tup[i],nodesets,&tup[dim+1],&node[dim+1]);CHKERRQ(ierr); 135 for (j = 0; j < dim+1; j++) node[j+(j>=i)] += wi * node[dim+1+j]; 136 w += wi; 137 } 138 for (i = 0; i < dim+1; i++) node[i] /= w; 139 } 140 PetscFunctionReturn(0); 141 } 142 143 /* compute simplex nodes for the biunit simplex from the 1D node family */ 144 static PetscErrorCode Petsc1DNodeFamilyComputeSimplexNodes(Petsc1DNodeFamily f, PetscInt dim, PetscInt degree, PetscReal points[]) 145 { 146 PetscInt *tup; 147 PetscInt k; 148 PetscInt npoints; 149 PetscReal **nodesets = NULL; 150 PetscInt worksize; 151 PetscReal *nodework; 152 PetscInt *tupwork; 153 PetscErrorCode ierr; 154 155 PetscFunctionBegin; 156 if (dim < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative dimension\n"); 157 if (degree < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have non-negative degree\n"); 158 if (!dim) PetscFunctionReturn(0); 159 ierr = PetscCalloc1(dim+2, &tup);CHKERRQ(ierr); 160 k = 0; 161 ierr = PetscDTBinomialInt(degree + dim, dim, &npoints);CHKERRQ(ierr); 162 ierr = Petsc1DNodeFamilyGetNodeSets(f, degree, &nodesets);CHKERRQ(ierr); 163 worksize = ((dim + 2) * (dim + 3)) / 2; 164 ierr = PetscMalloc2(worksize, &nodework, worksize, &tupwork);CHKERRQ(ierr); 165 for (k = 0; k < npoints; k++) { 166 PetscInt i; 167 168 tup[0] = degree; 169 for (i = 0; i < dim; i++) { 170 tup[0] -= tup[i+1]; 171 } 172 switch(f->nodeFamily) { 173 case PETSCDTNODES_EQUISPACED: 174 if (f->endpoints) { 175 for (i = 0; i < dim; i++) { 176 points[dim*k + i] = (PetscReal) tup[i+1] / (PetscReal) degree; 177 } 178 } else { 179 for (i = 0; i < dim; i++) { 180 points[dim*k + i] = ((PetscReal) tup[i+1] + 1./(dim+1.)) / (PetscReal) (degree + 1.); 181 } 182 } 183 break; 184 default: 185 for (i = 0; i < dim + 1; i++) tupwork[i] = tup[i]; 186 ierr = PetscNodeRecursive_Internal(dim, degree, nodesets, tupwork, nodework);CHKERRQ(ierr); 187 for (i = 0; i < dim; i++) points[dim*k + i] = nodework[i + 1]; 188 break; 189 } 190 ierr = PetscDualSpaceLatticePointLexicographic_Internal(dim, degree, &tup[1]);CHKERRQ(ierr); 191 } 192 /* map from unit simplex to biunit simplex */ 193 for (k = 0; k < npoints * dim; k++) points[k] = points[k] * 2. - 1.; 194 ierr = PetscFree2(nodework, tupwork);CHKERRQ(ierr); 195 ierr = PetscFree(tup); 196 PetscFunctionReturn(0); 197 } 198 199 struct _n_PetscLagNodeIndices 200 { 201 PetscInt refct; 202 PetscInt nodeIdxDim; 203 PetscInt nodeVecDim; 204 PetscInt nNodes; 205 PetscInt *nodeIdx; /* for each node an index of size nodeIdxDim */ 206 PetscReal *nodeVec; /* for each node a vector of size nodeVecDim */ 207 PetscInt *perm; /* if these are vertices, perm takes DMPlex point index to closure order; 208 if these are nodes, perm lists nodes in index revlex order */ 209 }; 210 211 PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices ni, PetscInt *nodeIdxDim, PetscInt *nodeVecDim, PetscInt *nNodes, const PetscInt *nodeIdx[], const PetscReal *nodeVec[]) 212 { 213 PetscFunctionBegin; 214 *nodeIdxDim = ni->nodeIdxDim; 215 *nodeVecDim = ni->nodeVecDim; 216 *nNodes = ni->nNodes; 217 *nodeIdx = ni->nodeIdx; 218 *nodeVec = ni->nodeVec; 219 PetscFunctionReturn(0); 220 } 221 222 static PetscErrorCode PetscLagNodeIndicesReference(PetscLagNodeIndices ni) 223 { 224 PetscFunctionBegin; 225 if (ni) ni->refct++; 226 PetscFunctionReturn(0); 227 } 228 229 static PetscErrorCode PetscLagNodeIndicesDestroy(PetscLagNodeIndices *ni) { 230 PetscErrorCode ierr; 231 232 PetscFunctionBegin; 233 if (!(*ni)) PetscFunctionReturn(0); 234 if (--(*ni)->refct > 0) { 235 *ni = NULL; 236 PetscFunctionReturn(0); 237 } 238 ierr = PetscFree((*ni)->nodeIdx);CHKERRQ(ierr); 239 ierr = PetscFree((*ni)->nodeVec);CHKERRQ(ierr); 240 ierr = PetscFree((*ni)->perm);CHKERRQ(ierr); 241 ierr = PetscFree(*ni);CHKERRQ(ierr); 242 *ni = NULL; 243 PetscFunctionReturn(0); 244 } 245 246 /* The vertex indices were written as though the vertices were in revlex order 247 * wrt coordinates. To understand the effect of different symmetries, we need 248 * them to be in closure order. We also need a permutation that takes point index 249 * to closure number */ 250 static PetscErrorCode PetscLagNodeIndicesComputeVertexOrder(DM dm, PetscLagNodeIndices ni, PetscBool sortIdx) 251 { 252 PetscInt v, w, vStart, vEnd, c, d; 253 PetscInt nVerts; 254 PetscInt closureSize = 0; 255 PetscInt *closure = NULL; 256 PetscInt *closureOrder; 257 PetscInt *invClosureOrder; 258 PetscInt *revlexOrder; 259 PetscInt *newNodeIdx; 260 PetscInt dim; 261 Vec coordVec; 262 const PetscScalar *coords; 263 PetscErrorCode ierr; 264 265 PetscFunctionBegin; 266 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 267 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 268 nVerts = vEnd - vStart; 269 ierr = PetscMalloc1(nVerts, &closureOrder);CHKERRQ(ierr); 270 ierr = PetscMalloc1(nVerts, &invClosureOrder);CHKERRQ(ierr); 271 ierr = PetscMalloc1(nVerts, &revlexOrder);CHKERRQ(ierr); 272 if (sortIdx) { 273 PetscInt nodeIdxDim = ni->nodeIdxDim; 274 PetscInt *idxOrder; 275 276 ierr = PetscMalloc1(nVerts * nodeIdxDim, &newNodeIdx);CHKERRQ(ierr); 277 ierr = PetscMalloc1(nVerts, &idxOrder);CHKERRQ(ierr); 278 for (v = 0; v < nVerts; v++) idxOrder[v] = v; 279 for (v = 0; v < nVerts; v++) { 280 for (w = v + 1; w < nVerts; w++) { 281 const PetscInt *iv = &(ni->nodeIdx[idxOrder[v] * nodeIdxDim]); 282 const PetscInt *iw = &(ni->nodeIdx[idxOrder[w] * nodeIdxDim]); 283 PetscInt diff = 0; 284 285 for (d = nodeIdxDim - 1; d >= 0; d--) if ((diff = (iv[d] - iw[d]))) break; 286 if (diff > 0) { 287 PetscInt swap = idxOrder[v]; 288 289 idxOrder[v] = idxOrder[w]; 290 idxOrder[w] = swap; 291 } 292 } 293 } 294 for (v = 0; v < nVerts; v++) { 295 for (d = 0; d < nodeIdxDim; d++) { 296 newNodeIdx[v * ni->nodeIdxDim + d] = ni->nodeIdx[idxOrder[v] * nodeIdxDim + d]; 297 } 298 } 299 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr); 300 ni->nodeIdx = newNodeIdx; 301 newNodeIdx = NULL; 302 ierr = PetscFree(idxOrder);CHKERRQ(ierr); 303 } 304 ierr = DMPlexGetTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 305 c = closureSize - nVerts; 306 for (v = 0; v < nVerts; v++) closureOrder[v] = closure[2 * (c + v)] - vStart; 307 for (v = 0; v < nVerts; v++) invClosureOrder[closureOrder[v]] = v; 308 ierr = DMPlexRestoreTransitiveClosure(dm, 0, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 309 ierr = DMGetCoordinatesLocal(dm, &coordVec);CHKERRQ(ierr); 310 ierr = VecGetArrayRead(coordVec, &coords);CHKERRQ(ierr); 311 /* bubble sort closure vertices by coordinates in revlex order */ 312 for (v = 0; v < nVerts; v++) revlexOrder[v] = v; 313 for (v = 0; v < nVerts; v++) { 314 for (w = v + 1; w < nVerts; w++) { 315 const PetscScalar *cv = &coords[closureOrder[revlexOrder[v]] * dim]; 316 const PetscScalar *cw = &coords[closureOrder[revlexOrder[w]] * dim]; 317 PetscReal diff = 0; 318 319 for (d = dim - 1; d >= 0; d--) if ((diff = PetscRealPart(cv[d] - cw[d])) != 0.) break; 320 if (diff > 0.) { 321 PetscInt swap = revlexOrder[v]; 322 323 revlexOrder[v] = revlexOrder[w]; 324 revlexOrder[w] = swap; 325 } 326 } 327 } 328 ierr = VecRestoreArrayRead(coordVec, &coords);CHKERRQ(ierr); 329 ierr = PetscMalloc1(ni->nodeIdxDim * nVerts, &newNodeIdx);CHKERRQ(ierr); 330 /* reorder nodeIdx to be in closure order */ 331 for (v = 0; v < nVerts; v++) { 332 for (d = 0; d < ni->nodeIdxDim; d++) { 333 newNodeIdx[revlexOrder[v] * ni->nodeIdxDim + d] = ni->nodeIdx[v * ni->nodeIdxDim + d]; 334 } 335 } 336 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr); 337 ni->nodeIdx = newNodeIdx; 338 ni->perm = invClosureOrder; 339 ierr = PetscFree(revlexOrder);CHKERRQ(ierr); 340 ierr = PetscFree(closureOrder);CHKERRQ(ierr); 341 PetscFunctionReturn(0); 342 } 343 344 static PetscErrorCode PetscLagNodeIndicesCreateSimplexVertices(DM dm, PetscLagNodeIndices *nodeIndices) 345 { 346 PetscLagNodeIndices ni; 347 PetscInt dim, d; 348 349 PetscErrorCode ierr; 350 351 PetscFunctionBegin; 352 ierr = PetscNew(&ni);CHKERRQ(ierr); 353 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 354 ni->nodeIdxDim = dim + 1; 355 ni->nodeVecDim = 0; 356 ni->nNodes = dim + 1; 357 ni->refct = 1; 358 ierr = PetscCalloc1((dim + 1)*(dim + 1), &(ni->nodeIdx));CHKERRQ(ierr); 359 for (d = 0; d < dim + 1; d++) ni->nodeIdx[d*(dim + 2)] = 1; 360 ierr = PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_FALSE);CHKERRQ(ierr); 361 *nodeIndices = ni; 362 PetscFunctionReturn(0); 363 } 364 365 static PetscErrorCode PetscLagNodeIndicesCreateTensorVertices(DM dm, PetscLagNodeIndices facetni, PetscLagNodeIndices *nodeIndices) 366 { 367 PetscLagNodeIndices ni; 368 PetscInt nodeIdxDim, subNodeIdxDim = facetni->nodeIdxDim; 369 PetscInt nVerts, nSubVerts = facetni->nNodes; 370 PetscInt dim, d, e, f, g; 371 372 PetscErrorCode ierr; 373 374 PetscFunctionBegin; 375 ierr = PetscNew(&ni);CHKERRQ(ierr); 376 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 377 ni->nodeIdxDim = nodeIdxDim = subNodeIdxDim + 2; 378 ni->nodeVecDim = 0; 379 ni->nNodes = nVerts = 2 * nSubVerts; 380 ni->refct = 1; 381 ierr = PetscCalloc1(nodeIdxDim * nVerts, &(ni->nodeIdx));CHKERRQ(ierr); 382 for (f = 0, d = 0; d < 2; d++) { 383 for (e = 0; e < nSubVerts; e++, f++) { 384 for (g = 0; g < subNodeIdxDim; g++) { 385 ni->nodeIdx[f * nodeIdxDim + g] = facetni->nodeIdx[e * subNodeIdxDim + g]; 386 } 387 ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim] = (1 - d); 388 ni->nodeIdx[f * nodeIdxDim + subNodeIdxDim + 1] = d; 389 } 390 } 391 ierr = PetscLagNodeIndicesComputeVertexOrder(dm, ni, PETSC_TRUE);CHKERRQ(ierr); 392 *nodeIndices = ni; 393 PetscFunctionReturn(0); 394 } 395 396 static PetscErrorCode PetscLagNodeIndicesPushForward(DM dm, PetscLagNodeIndices vert, PetscInt p, PetscLagNodeIndices vertp, PetscLagNodeIndices nodep, PetscInt ornt, PetscInt formDegree, PetscInt pfNodeIdx[], PetscReal pfNodeVec[]) 397 { 398 PetscInt *closureVerts; 399 PetscInt closureSize = 0; 400 PetscInt *closure = NULL; 401 PetscInt dim, pdim, c, i, j, k, n, v, vStart, vEnd; 402 PetscInt nSubVert = vertp->nNodes; 403 PetscInt nodeIdxDim = vert->nodeIdxDim; 404 PetscInt subNodeIdxDim = vertp->nodeIdxDim; 405 PetscInt nNodes = nodep->nNodes; 406 const PetscInt *vertIdx = vert->nodeIdx; 407 const PetscInt *subVertIdx = vertp->nodeIdx; 408 const PetscInt *nodeIdx = nodep->nodeIdx; 409 const PetscReal *nodeVec = nodep->nodeVec; 410 PetscReal *J, *Jstar; 411 PetscReal detJ; 412 PetscInt depth, pdepth, Nk, pNk; 413 Vec coordVec; 414 PetscScalar *newCoords = NULL; 415 const PetscScalar *oldCoords = NULL; 416 PetscErrorCode ierr; 417 418 PetscFunctionBegin; 419 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 420 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 421 ierr = DMGetCoordinatesLocal(dm, &coordVec);CHKERRQ(ierr); 422 ierr = DMPlexGetPointDepth(dm, p, &pdepth);CHKERRQ(ierr); 423 pdim = pdepth != depth ? pdepth != 0 ? pdepth : 0 : dim; 424 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 425 ierr = DMGetWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);CHKERRQ(ierr); 426 ierr = DMPlexGetTransitiveClosure_Internal(dm, p, ornt, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 427 c = closureSize - nSubVert; 428 /* we want which cell closure indices the closure of this point corresponds to */ 429 for (v = 0; v < nSubVert; v++) closureVerts[v] = vert->perm[closure[2 * (c + v)] - vStart]; 430 ierr = DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 431 /* push forward indices */ 432 for (i = 0; i < nodeIdxDim; i++) { /* for every component of the target index space */ 433 /* check if this is a component that all vertices around this point have in common */ 434 for (j = 1; j < nSubVert; j++) { 435 if (vertIdx[closureVerts[j] * nodeIdxDim + i] != vertIdx[closureVerts[0] * nodeIdxDim + i]) break; 436 } 437 if (j == nSubVert) { /* all vertices have this component in common, directly copy to output */ 438 PetscInt val = vertIdx[closureVerts[0] * nodeIdxDim + i]; 439 for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = val; 440 } else { 441 PetscInt subi = -1; 442 /* there must be a component in vertp that looks the same */ 443 for (k = 0; k < subNodeIdxDim; k++) { 444 for (j = 0; j < nSubVert; j++) { 445 if (vertIdx[closureVerts[j] * nodeIdxDim + i] != subVertIdx[j * subNodeIdxDim + k]) break; 446 } 447 if (j == nSubVert) { 448 subi = k; 449 break; 450 } 451 } 452 if (subi < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Did not find matching coordinate\n"); 453 for (n = 0; n < nNodes; n++) pfNodeIdx[n * nodeIdxDim + i] = nodeIdx[n * subNodeIdxDim + subi]; 454 } 455 } 456 /* push forward vectors */ 457 ierr = DMGetWorkArray(dm, dim * dim, MPIU_REAL, &J);CHKERRQ(ierr); 458 if (ornt != 0) { 459 PetscInt closureSize2 = 0; 460 PetscInt *closure2 = NULL; 461 462 ierr = DMPlexGetTransitiveClosure_Internal(dm, p, 0, PETSC_TRUE, &closureSize2, &closure2);CHKERRQ(ierr); 463 ierr = PetscMalloc1(dim * nSubVert, &newCoords);CHKERRQ(ierr); 464 ierr = VecGetArrayRead(coordVec, &oldCoords);CHKERRQ(ierr); 465 for (v = 0; v < nSubVert; v++) { 466 PetscInt d; 467 for (d = 0; d < dim; d++) { 468 newCoords[(closure2[2 * (c + v)] - vStart) * dim + d] = oldCoords[closureVerts[v] * dim + d]; 469 } 470 } 471 ierr = VecRestoreArrayRead(coordVec, &oldCoords);CHKERRQ(ierr); 472 ierr = DMPlexRestoreTransitiveClosure(dm, p, PETSC_TRUE, &closureSize2, &closure2);CHKERRQ(ierr); 473 ierr = VecPlaceArray(coordVec, newCoords);CHKERRQ(ierr); 474 } 475 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, NULL, J, NULL, &detJ);CHKERRQ(ierr); 476 if (ornt != 0) { 477 ierr = VecResetArray(coordVec);CHKERRQ(ierr); 478 ierr = PetscFree(newCoords);CHKERRQ(ierr); 479 } 480 ierr = DMRestoreWorkArray(dm, nSubVert, MPIU_INT, &closureVerts);CHKERRQ(ierr); 481 /* compactify */ 482 for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j]; 483 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr); 484 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(formDegree), &pNk);CHKERRQ(ierr); 485 ierr = DMGetWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);CHKERRQ(ierr); 486 ierr = PetscDTAltVPullbackMatrix(pdim, dim, J, formDegree, Jstar);CHKERRQ(ierr); 487 for (n = 0; n < nNodes; n++) { 488 for (i = 0; i < Nk; i++) { 489 PetscReal val = 0.; 490 for (j = 0; j < pNk; j++) val += nodeVec[n * pNk + j] * Jstar[j * pNk + i]; 491 pfNodeVec[n * Nk + i] = val; 492 } 493 } 494 ierr = DMRestoreWorkArray(dm, pNk * Nk, MPIU_REAL, &Jstar);CHKERRQ(ierr); 495 ierr = DMRestoreWorkArray(dm, dim * dim, MPIU_REAL, &J);CHKERRQ(ierr); 496 PetscFunctionReturn(0); 497 } 498 499 static PetscErrorCode PetscLagNodeIndicesTensor(PetscLagNodeIndices tracei, PetscInt dimT, PetscInt kT, PetscLagNodeIndices fiberi, PetscInt dimF, PetscInt kF, PetscLagNodeIndices *nodeIndices) 500 { 501 PetscInt dim = dimT + dimF; 502 PetscInt nodeIdxDim, nNodes; 503 PetscInt formDegree = kT + kF; 504 PetscInt Nk, NkT, NkF; 505 PetscInt MkT, MkF; 506 PetscLagNodeIndices ni; 507 PetscInt i, j, l; 508 PetscReal *projF, *projT; 509 PetscReal *projFstar, *projTstar; 510 PetscReal *workF, *workF2, *workT, *workT2, *work, *work2; 511 PetscReal *wedgeMat; 512 PetscReal sign; 513 PetscErrorCode ierr; 514 515 PetscFunctionBegin; 516 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr); 517 ierr = PetscDTBinomialInt(dimT, PetscAbsInt(kT), &NkT);CHKERRQ(ierr); 518 ierr = PetscDTBinomialInt(dimF, PetscAbsInt(kF), &NkF);CHKERRQ(ierr); 519 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kT), &MkT);CHKERRQ(ierr); 520 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kF), &MkF);CHKERRQ(ierr); 521 ierr = PetscNew(&ni);CHKERRQ(ierr); 522 ni->nodeIdxDim = nodeIdxDim = tracei->nodeIdxDim + fiberi->nodeIdxDim; 523 ni->nodeVecDim = Nk; 524 ni->nNodes = nNodes = tracei->nNodes * fiberi->nNodes; 525 ni->refct = 1; 526 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr); 527 /* first concatenate the indices */ 528 for (l = 0, j = 0; j < fiberi->nNodes; j++) { 529 for (i = 0; i < tracei->nNodes; i++, l++) { 530 PetscInt m, n = 0; 531 532 for (m = 0; m < tracei->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = tracei->nodeIdx[i * tracei->nodeIdxDim + m]; 533 for (m = 0; m < fiberi->nodeIdxDim; m++) ni->nodeIdx[l * nodeIdxDim + n++] = fiberi->nodeIdx[j * fiberi->nodeIdxDim + m]; 534 } 535 } 536 537 /* now wedge together the push-forward vectors */ 538 ierr = PetscMalloc1(nNodes * Nk, &(ni->nodeVec));CHKERRQ(ierr); 539 ierr = PetscCalloc2(dimT*dim, &projT, dimF*dim, &projF);CHKERRQ(ierr); 540 for (i = 0; i < dimT; i++) projT[i * (dim + 1)] = 1.; 541 for (i = 0; i < dimF; i++) projF[i * (dim + dimT + 1) + dimT] = 1.; 542 ierr = PetscMalloc2(MkT*NkT, &projTstar, MkF*NkF, &projFstar);CHKERRQ(ierr); 543 ierr = PetscDTAltVPullbackMatrix(dim, dimT, projT, kT, projTstar);CHKERRQ(ierr); 544 ierr = PetscDTAltVPullbackMatrix(dim, dimF, projF, kF, projFstar);CHKERRQ(ierr); 545 ierr = PetscMalloc6(MkT, &workT, MkT, &workT2, MkF, &workF, MkF, &workF2, Nk, &work, Nk, &work2);CHKERRQ(ierr); 546 ierr = PetscMalloc1(Nk * MkT, &wedgeMat);CHKERRQ(ierr); 547 sign = (PetscAbsInt(kT * kF) & 1) ? -1. : 1.; 548 for (l = 0, j = 0; j < fiberi->nNodes; j++) { 549 PetscInt d, e; 550 551 /* push forward fiber k-form */ 552 for (d = 0; d < MkF; d++) { 553 PetscReal val = 0.; 554 for (e = 0; e < NkF; e++) val += projFstar[d * NkF + e] * fiberi->nodeVec[j * NkF + e]; 555 workF[d] = val; 556 } 557 /* Hodge star to proper form if necessary */ 558 if (kF < 0) { 559 for (d = 0; d < MkF; d++) workF2[d] = workF[d]; 560 ierr = PetscDTAltVStar(dim, PetscAbsInt(kF), 1, workF2, workF);CHKERRQ(ierr); 561 } 562 /* Compute the matrix that wedges this form with one of the trace k-form */ 563 ierr = PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kF), PetscAbsInt(kT), workF, wedgeMat);CHKERRQ(ierr); 564 for (i = 0; i < tracei->nNodes; i++, l++) { 565 /* push forward trace k-form */ 566 for (d = 0; d < MkT; d++) { 567 PetscReal val = 0.; 568 for (e = 0; e < NkT; e++) val += projTstar[d * NkT + e] * tracei->nodeVec[i * NkT + e]; 569 workT[d] = val; 570 } 571 /* Hodge star to proper form if necessary */ 572 if (kT < 0) { 573 for (d = 0; d < MkT; d++) workT2[d] = workT[d]; 574 ierr = PetscDTAltVStar(dim, PetscAbsInt(kT), 1, workT2, workT);CHKERRQ(ierr); 575 } 576 /* compute the wedge product of the push-forward trace form and firer forms */ 577 for (d = 0; d < Nk; d++) { 578 PetscReal val = 0.; 579 for (e = 0; e < MkT; e++) val += wedgeMat[d * MkT + e] * workT[e]; 580 work[d] = val; 581 } 582 /* inverse Hodge star from proper form if necessary */ 583 if (formDegree < 0) { 584 for (d = 0; d < Nk; d++) work2[d] = work[d]; 585 ierr = PetscDTAltVStar(dim, PetscAbsInt(formDegree), -1, work2, work);CHKERRQ(ierr); 586 } 587 /* insert into the array (adjusting for sign) */ 588 for (d = 0; d < Nk; d++) ni->nodeVec[l * Nk + d] = sign * work[d]; 589 } 590 } 591 ierr = PetscFree(wedgeMat);CHKERRQ(ierr); 592 ierr = PetscFree6(workT, workT2, workF, workF2, work, work2);CHKERRQ(ierr); 593 ierr = PetscFree2(projTstar, projFstar);CHKERRQ(ierr); 594 ierr = PetscFree2(projT, projF);CHKERRQ(ierr); 595 *nodeIndices = ni; 596 PetscFunctionReturn(0); 597 } 598 599 static PetscErrorCode PetscLagNodeIndicesMerge(PetscLagNodeIndices niA, PetscLagNodeIndices niB, PetscLagNodeIndices *nodeIndices) 600 { 601 PetscLagNodeIndices ni; 602 PetscInt nodeIdxDim, nodeVecDim, nNodes; 603 PetscErrorCode ierr; 604 605 PetscFunctionBegin; 606 ierr = PetscNew(&ni);CHKERRQ(ierr); 607 ni->nodeIdxDim = nodeIdxDim = niA->nodeIdxDim; 608 if (niB->nodeIdxDim != nodeIdxDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeIdxDim"); 609 ni->nodeVecDim = nodeVecDim = niA->nodeVecDim; 610 if (niB->nodeVecDim != nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Cannot merge PetscLagNodeIndices with different nodeVecDim"); 611 ni->nNodes = nNodes = niA->nNodes + niB->nNodes; 612 ni->refct = 1; 613 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr); 614 ierr = PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));CHKERRQ(ierr); 615 ierr = PetscArraycpy(ni->nodeIdx, niA->nodeIdx, niA->nNodes * nodeIdxDim);CHKERRQ(ierr); 616 ierr = PetscArraycpy(ni->nodeVec, niA->nodeVec, niA->nNodes * nodeVecDim);CHKERRQ(ierr); 617 ierr = PetscArraycpy(&(ni->nodeIdx[niA->nNodes * nodeIdxDim]), niB->nodeIdx, niB->nNodes * nodeIdxDim);CHKERRQ(ierr); 618 ierr = PetscArraycpy(&(ni->nodeVec[niA->nNodes * nodeVecDim]), niB->nodeVec, niB->nNodes * nodeVecDim);CHKERRQ(ierr); 619 *nodeIndices = ni; 620 PetscFunctionReturn(0); 621 } 622 623 #define PETSCTUPINTCOMPREVLEX(N) \ 624 static int PetscTupIntCompRevlex_##N(const void *a, const void *b) \ 625 { \ 626 const PetscInt *A = (const PetscInt *) a; \ 627 const PetscInt *B = (const PetscInt *) b; \ 628 int i; \ 629 PetscInt diff = 0; \ 630 for (i = 0; i < N; i++) { \ 631 diff = A[N - i] - B[N - i]; \ 632 if (diff) break; \ 633 } \ 634 return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1; \ 635 } 636 637 PETSCTUPINTCOMPREVLEX(3) 638 PETSCTUPINTCOMPREVLEX(4) 639 PETSCTUPINTCOMPREVLEX(5) 640 PETSCTUPINTCOMPREVLEX(6) 641 PETSCTUPINTCOMPREVLEX(7) 642 643 static int PetscTupIntCompRevlex_N(const void *a, const void *b) 644 { 645 const PetscInt *A = (const PetscInt *) a; 646 const PetscInt *B = (const PetscInt *) b; 647 int i; 648 int N = A[0]; 649 PetscInt diff = 0; 650 for (i = 0; i < N; i++) { 651 diff = A[N - i] - B[N - i]; 652 if (diff) break; 653 } 654 return (diff <= 0) ? (diff < 0) ? -1 : 0 : 1; 655 } 656 657 static PetscErrorCode PetscLagNodeIndicesGetPermutation(PetscLagNodeIndices ni, PetscInt *perm[]) 658 { 659 PetscErrorCode ierr; 660 661 PetscFunctionBegin; 662 if (!(ni->perm)) { 663 PetscInt *sorter; 664 PetscInt m = ni->nNodes; 665 PetscInt nodeIdxDim = ni->nodeIdxDim; 666 PetscInt i, j, k, l; 667 PetscInt *prm; 668 int (*comp) (const void *, const void *); 669 670 ierr = PetscMalloc1((nodeIdxDim + 2) * m, &sorter);CHKERRQ(ierr); 671 for (k = 0, l = 0, i = 0; i < m; i++) { 672 sorter[k++] = nodeIdxDim + 1; 673 sorter[k++] = i; 674 for (j = 0; j < nodeIdxDim; j++) sorter[k++] = ni->nodeIdx[l++]; 675 } 676 switch (nodeIdxDim) { 677 case 2: 678 comp = PetscTupIntCompRevlex_3; 679 break; 680 case 3: 681 comp = PetscTupIntCompRevlex_4; 682 break; 683 case 4: 684 comp = PetscTupIntCompRevlex_5; 685 break; 686 case 5: 687 comp = PetscTupIntCompRevlex_6; 688 break; 689 case 6: 690 comp = PetscTupIntCompRevlex_7; 691 break; 692 default: 693 comp = PetscTupIntCompRevlex_N; 694 break; 695 } 696 qsort(sorter, m, (nodeIdxDim + 2) * sizeof(PetscInt), comp); 697 ierr = PetscMalloc1(m, &prm);CHKERRQ(ierr); 698 for (i = 0; i < m; i++) prm[i] = sorter[(nodeIdxDim + 2) * i + 1]; 699 ni->perm = prm; 700 ierr = PetscFree(sorter); 701 } 702 *perm = ni->perm; 703 PetscFunctionReturn(0); 704 } 705 706 static PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp) 707 { 708 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 709 PetscErrorCode ierr; 710 711 PetscFunctionBegin; 712 if (lag->symperms) { 713 PetscInt **selfSyms = lag->symperms[0]; 714 715 if (selfSyms) { 716 PetscInt i, **allocated = &selfSyms[-lag->selfSymOff]; 717 718 for (i = 0; i < lag->numSelfSym; i++) { 719 ierr = PetscFree(allocated[i]);CHKERRQ(ierr); 720 } 721 ierr = PetscFree(allocated);CHKERRQ(ierr); 722 } 723 ierr = PetscFree(lag->symperms);CHKERRQ(ierr); 724 } 725 if (lag->symflips) { 726 PetscScalar **selfSyms = lag->symflips[0]; 727 728 if (selfSyms) { 729 PetscInt i; 730 PetscScalar **allocated = &selfSyms[-lag->selfSymOff]; 731 732 for (i = 0; i < lag->numSelfSym; i++) { 733 ierr = PetscFree(allocated[i]);CHKERRQ(ierr); 734 } 735 ierr = PetscFree(allocated);CHKERRQ(ierr); 736 } 737 ierr = PetscFree(lag->symflips);CHKERRQ(ierr); 738 } 739 ierr = Petsc1DNodeFamilyDestroy(&(lag->nodeFamily));CHKERRQ(ierr); 740 ierr = PetscLagNodeIndicesDestroy(&(lag->vertIndices));CHKERRQ(ierr); 741 ierr = PetscLagNodeIndicesDestroy(&(lag->intNodeIndices));CHKERRQ(ierr); 742 ierr = PetscLagNodeIndicesDestroy(&(lag->allNodeIndices));CHKERRQ(ierr); 743 ierr = PetscFree(lag);CHKERRQ(ierr); 744 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);CHKERRQ(ierr); 745 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);CHKERRQ(ierr); 746 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", NULL);CHKERRQ(ierr); 747 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", NULL);CHKERRQ(ierr); 748 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL);CHKERRQ(ierr); 749 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL);CHKERRQ(ierr); 750 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL);CHKERRQ(ierr); 751 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL);CHKERRQ(ierr); 752 PetscFunctionReturn(0); 753 } 754 755 static PetscErrorCode PetscDualSpaceLagrangeView_Ascii(PetscDualSpace sp, PetscViewer viewer) 756 { 757 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 758 PetscErrorCode ierr; 759 760 PetscFunctionBegin; 761 ierr = PetscViewerASCIIPrintf(viewer, "%s %s%sLagrange dual space\n", lag->continuous ? "Continuous" : "Discontinuous", lag->tensorSpace ? "tensor " : "", lag->trimmed ? "trimmed " : "");CHKERRQ(ierr); 762 PetscFunctionReturn(0); 763 } 764 765 static PetscErrorCode PetscDualSpaceView_Lagrange(PetscDualSpace sp, PetscViewer viewer) 766 { 767 PetscBool iascii; 768 PetscErrorCode ierr; 769 770 PetscFunctionBegin; 771 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 772 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2); 773 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 774 if (iascii) {ierr = PetscDualSpaceLagrangeView_Ascii(sp, viewer);CHKERRQ(ierr);} 775 PetscFunctionReturn(0); 776 } 777 778 static PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp) 779 { 780 PetscBool continuous, tensor, trimmed, flg, flg2, flg3; 781 PetscDTNodeType nodeType; 782 PetscReal nodeExponent; 783 PetscBool nodeEndpoints; 784 PetscErrorCode ierr; 785 786 PetscFunctionBegin; 787 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &continuous);CHKERRQ(ierr); 788 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr); 789 ierr = PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);CHKERRQ(ierr); 790 ierr = PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &nodeEndpoints, &nodeExponent);CHKERRQ(ierr); 791 if (nodeType == PETSCDTNODES_DEFAULT) nodeType = PETSCDTNODES_GAUSSJACOBI; 792 ierr = PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");CHKERRQ(ierr); 793 ierr = PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);CHKERRQ(ierr); 794 if (flg) {ierr = PetscDualSpaceLagrangeSetContinuity(sp, continuous);CHKERRQ(ierr);} 795 ierr = PetscOptionsBool("-petscdualspace_lagrange_tensor", "Flag for tensor dual space", "PetscDualSpaceLagrangeSetTensor", tensor, &tensor, &flg);CHKERRQ(ierr); 796 if (flg) {ierr = PetscDualSpaceLagrangeSetTensor(sp, tensor);CHKERRQ(ierr);} 797 ierr = PetscOptionsBool("-petscdualspace_lagrange_trimmed", "Flag for trimmed dual space", "PetscDualSpaceLagrangeSetTrimmed", trimmed, &trimmed, &flg);CHKERRQ(ierr); 798 if (flg) {ierr = PetscDualSpaceLagrangeSetTrimmed(sp, trimmed);CHKERRQ(ierr);} 799 ierr = PetscOptionsEnum("-petscdualspace_lagrange_node_type", "Lagrange node location type", "PetscDualSpaceLagrangeSetNodeType", PetscDTNodeTypes, (PetscEnum)nodeType, (PetscEnum *)&nodeType, &flg);CHKERRQ(ierr); 800 ierr = PetscOptionsBool("-petscdualspace_lagrange_node_endpoints", "Flag for nodes that include endpoints", "PetscDualSpaceLagrangeSetNodeType", nodeEndpoints, &nodeEndpoints, &flg2);CHKERRQ(ierr); 801 flg3 = PETSC_FALSE; 802 if (nodeType == PETSCDTNODES_GAUSSJACOBI) { 803 ierr = PetscOptionsReal("-petscdualspace_lagrange_node_exponent", "Gauss-Jacobi weight function exponent", "PetscDualSpaceLagrangeSetNodeType", nodeExponent, &nodeExponent, &flg3);CHKERRQ(ierr); 804 } 805 if (flg || flg2 || flg3) {ierr = PetscDualSpaceLagrangeSetNodeType(sp, nodeType, nodeEndpoints, nodeExponent);CHKERRQ(ierr);} 806 ierr = PetscOptionsTail();CHKERRQ(ierr); 807 PetscFunctionReturn(0); 808 } 809 810 static PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace spNew) 811 { 812 PetscBool cont, tensor, trimmed, boundary; 813 PetscDTNodeType nodeType; 814 PetscReal exponent; 815 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 816 PetscErrorCode ierr; 817 818 PetscFunctionBegin; 819 ierr = PetscDualSpaceLagrangeGetContinuity(sp, &cont);CHKERRQ(ierr); 820 ierr = PetscDualSpaceLagrangeSetContinuity(spNew, cont);CHKERRQ(ierr); 821 ierr = PetscDualSpaceLagrangeGetTensor(sp, &tensor);CHKERRQ(ierr); 822 ierr = PetscDualSpaceLagrangeSetTensor(spNew, tensor);CHKERRQ(ierr); 823 ierr = PetscDualSpaceLagrangeGetTrimmed(sp, &trimmed);CHKERRQ(ierr); 824 ierr = PetscDualSpaceLagrangeSetTrimmed(spNew, trimmed);CHKERRQ(ierr); 825 ierr = PetscDualSpaceLagrangeGetNodeType(sp, &nodeType, &boundary, &exponent);CHKERRQ(ierr); 826 ierr = PetscDualSpaceLagrangeSetNodeType(spNew, nodeType, boundary, exponent);CHKERRQ(ierr); 827 if (lag->nodeFamily) { 828 PetscDualSpace_Lag *lagnew = (PetscDualSpace_Lag *) spNew->data; 829 830 ierr = Petsc1DNodeFamilyReference(lag->nodeFamily);CHKERRQ(ierr); 831 lagnew->nodeFamily = lag->nodeFamily; 832 } 833 PetscFunctionReturn(0); 834 } 835 836 static PetscErrorCode PetscDualSpaceCreateEdgeSubspace_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt k, PetscInt Nc, PetscBool interiorOnly, PetscDualSpace *bdsp) 837 { 838 DM K; 839 PetscDualSpace_Lag *newlag; 840 PetscErrorCode ierr; 841 842 PetscFunctionBegin; 843 ierr = PetscDualSpaceDuplicate(sp,bdsp);CHKERRQ(ierr); 844 ierr = PetscDualSpaceSetFormDegree(*bdsp, k);CHKERRQ(ierr); 845 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, 1, PETSC_TRUE, &K);CHKERRQ(ierr); 846 ierr = PetscDualSpaceSetDM(*bdsp, K);CHKERRQ(ierr); 847 ierr = DMDestroy(&K);CHKERRQ(ierr); 848 ierr = PetscDualSpaceSetOrder(*bdsp, order);CHKERRQ(ierr); 849 ierr = PetscDualSpaceSetNumComponents(*bdsp, Nc);CHKERRQ(ierr); 850 newlag = (PetscDualSpace_Lag *) (*bdsp)->data; 851 newlag->interiorOnly = interiorOnly; 852 ierr = PetscDualSpaceSetUp(*bdsp);CHKERRQ(ierr); 853 PetscFunctionReturn(0); 854 } 855 856 /* just the points, weights aren't handled */ 857 static PetscErrorCode PetscQuadratureCreateTensor(PetscQuadrature trace, PetscQuadrature fiber, PetscQuadrature *product) 858 { 859 PetscInt dimTrace, dimFiber; 860 PetscInt numPointsTrace, numPointsFiber; 861 PetscInt dim, numPoints; 862 const PetscReal *pointsTrace; 863 const PetscReal *pointsFiber; 864 PetscReal *points; 865 PetscInt i, j, k, p; 866 PetscErrorCode ierr; 867 868 PetscFunctionBegin; 869 ierr = PetscQuadratureGetData(trace, &dimTrace, NULL, &numPointsTrace, &pointsTrace, NULL);CHKERRQ(ierr); 870 ierr = PetscQuadratureGetData(fiber, &dimFiber, NULL, &numPointsFiber, &pointsFiber, NULL);CHKERRQ(ierr); 871 dim = dimTrace + dimFiber; 872 numPoints = numPointsFiber * numPointsTrace; 873 ierr = PetscMalloc1(numPoints * dim, &points);CHKERRQ(ierr); 874 for (p = 0, j = 0; j < numPointsFiber; j++) { 875 for (i = 0; i < numPointsTrace; i++, p++) { 876 for (k = 0; k < dimTrace; k++) points[p * dim + k] = pointsTrace[i * dimTrace + k]; 877 for (k = 0; k < dimFiber; k++) points[p * dim + dimTrace + k] = pointsFiber[j * dimFiber + k]; 878 } 879 } 880 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, product);CHKERRQ(ierr); 881 ierr = PetscQuadratureSetData(*product, dim, 0, numPoints, points, NULL);CHKERRQ(ierr); 882 PetscFunctionReturn(0); 883 } 884 885 static PetscErrorCode MatTensorAltV(Mat trace, Mat fiber, PetscInt dimTrace, PetscInt kTrace, PetscInt dimFiber, PetscInt kFiber, Mat *product) 886 { 887 PetscInt mTrace, nTrace, mFiber, nFiber, m, n, k, i, j, l; 888 PetscInt dim, NkTrace, NkFiber, Nk; 889 PetscInt dT, dF; 890 PetscInt *nnzTrace, *nnzFiber, *nnz; 891 PetscInt iT, iF, jT, jF, il, jl; 892 PetscReal *workT, *workT2, *workF, *workF2, *work, *workstar; 893 PetscReal *projT, *projF; 894 PetscReal *projTstar, *projFstar; 895 PetscReal *wedgeMat; 896 PetscReal sign; 897 PetscScalar *workS; 898 Mat prod; 899 /* this produces dof groups that look like the identity */ 900 PetscErrorCode ierr; 901 902 PetscFunctionBegin; 903 ierr = MatGetSize(trace, &mTrace, &nTrace);CHKERRQ(ierr); 904 ierr = PetscDTBinomialInt(dimTrace, PetscAbsInt(kTrace), &NkTrace);CHKERRQ(ierr); 905 if (nTrace % NkTrace) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of trace matrix is not a multiple of k-form size"); 906 ierr = MatGetSize(fiber, &mFiber, &nFiber);CHKERRQ(ierr); 907 ierr = PetscDTBinomialInt(dimFiber, PetscAbsInt(kFiber), &NkFiber);CHKERRQ(ierr); 908 if (nFiber % NkFiber) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "point value space of fiber matrix is not a multiple of k-form size"); 909 ierr = PetscMalloc2(mTrace, &nnzTrace, mFiber, &nnzFiber);CHKERRQ(ierr); 910 for (i = 0; i < mTrace; i++) { 911 ierr = MatGetRow(trace, i, &(nnzTrace[i]), NULL, NULL);CHKERRQ(ierr); 912 if (nnzTrace[i] % NkTrace) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in trace matrix are not in k-form size blocks"); 913 } 914 for (i = 0; i < mFiber; i++) { 915 ierr = MatGetRow(fiber, i, &(nnzFiber[i]), NULL, NULL);CHKERRQ(ierr); 916 if (nnzFiber[i] % NkFiber) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in fiber matrix are not in k-form size blocks"); 917 } 918 dim = dimTrace + dimFiber; 919 k = kFiber + kTrace; 920 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 921 m = mTrace * mFiber; 922 ierr = PetscMalloc1(m, &nnz);CHKERRQ(ierr); 923 for (l = 0, j = 0; j < mFiber; j++) for (i = 0; i < mTrace; i++, l++) nnz[l] = (nnzTrace[i] / NkTrace) * (nnzFiber[j] / NkFiber) * Nk; 924 n = (nTrace / NkTrace) * (nFiber / NkFiber) * Nk; 925 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &prod);CHKERRQ(ierr); 926 ierr = PetscFree(nnz);CHKERRQ(ierr); 927 ierr = PetscFree2(nnzTrace,nnzFiber);CHKERRQ(ierr); 928 /* reasoning about which points each dof needs depends on having zeros computed at points preserved */ 929 ierr = MatSetOption(prod, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr); 930 /* compute pullbacks */ 931 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kTrace), &dT);CHKERRQ(ierr); 932 ierr = PetscDTBinomialInt(dim, PetscAbsInt(kFiber), &dF);CHKERRQ(ierr); 933 ierr = PetscMalloc4(dimTrace * dim, &projT, dimFiber * dim, &projF, dT * NkTrace, &projTstar, dF * NkFiber, &projFstar);CHKERRQ(ierr); 934 ierr = PetscArrayzero(projT, dimTrace * dim);CHKERRQ(ierr); 935 for (i = 0; i < dimTrace; i++) projT[i * (dim + 1)] = 1.; 936 ierr = PetscArrayzero(projF, dimFiber * dim);CHKERRQ(ierr); 937 for (i = 0; i < dimFiber; i++) projF[i * (dim + 1) + dimTrace] = 1.; 938 ierr = PetscDTAltVPullbackMatrix(dim, dimTrace, projT, kTrace, projTstar);CHKERRQ(ierr); 939 ierr = PetscDTAltVPullbackMatrix(dim, dimFiber, projF, kFiber, projFstar);CHKERRQ(ierr); 940 ierr = PetscMalloc5(dT, &workT, dF, &workF, Nk, &work, Nk, &workstar, Nk, &workS);CHKERRQ(ierr); 941 ierr = PetscMalloc2(dT, &workT2, dF, &workF2);CHKERRQ(ierr); 942 ierr = PetscMalloc1(Nk * dT, &wedgeMat);CHKERRQ(ierr); 943 sign = (PetscAbsInt(kTrace * kFiber) & 1) ? -1. : 1.; 944 for (i = 0, iF = 0; iF < mFiber; iF++) { 945 PetscInt ncolsF, nformsF; 946 const PetscInt *colsF; 947 const PetscScalar *valsF; 948 949 ierr = MatGetRow(fiber, iF, &ncolsF, &colsF, &valsF);CHKERRQ(ierr); 950 nformsF = ncolsF / NkFiber; 951 for (iT = 0; iT < mTrace; iT++, i++) { 952 PetscInt ncolsT, nformsT; 953 const PetscInt *colsT; 954 const PetscScalar *valsT; 955 956 ierr = MatGetRow(trace, iT, &ncolsT, &colsT, &valsT);CHKERRQ(ierr); 957 nformsT = ncolsT / NkTrace; 958 for (j = 0, jF = 0; jF < nformsF; jF++) { 959 PetscInt colF = colsF[jF * NkFiber] / NkFiber; 960 961 for (il = 0; il < dF; il++) { 962 PetscReal val = 0.; 963 for (jl = 0; jl < NkFiber; jl++) val += projFstar[il * NkFiber + jl] * PetscRealPart(valsF[jF * NkFiber + jl]); 964 workF[il] = val; 965 } 966 if (kFiber < 0) { 967 for (il = 0; il < dF; il++) workF2[il] = workF[il]; 968 ierr = PetscDTAltVStar(dim, PetscAbsInt(kFiber), 1, workF2, workF);CHKERRQ(ierr); 969 } 970 ierr = PetscDTAltVWedgeMatrix(dim, PetscAbsInt(kFiber), PetscAbsInt(kTrace), workF, wedgeMat);CHKERRQ(ierr); 971 for (jT = 0; jT < nformsT; jT++, j++) { 972 PetscInt colT = colsT[jT * NkTrace] / NkTrace; 973 PetscInt col = colF * (nTrace / NkTrace) + colT; 974 const PetscScalar *vals; 975 976 for (il = 0; il < dT; il++) { 977 PetscReal val = 0.; 978 for (jl = 0; jl < NkTrace; jl++) val += projTstar[il * NkTrace + jl] * PetscRealPart(valsT[jT * NkTrace + jl]); 979 workT[il] = val; 980 } 981 if (kTrace < 0) { 982 for (il = 0; il < dT; il++) workT2[il] = workT[il]; 983 ierr = PetscDTAltVStar(dim, PetscAbsInt(kTrace), 1, workT2, workT);CHKERRQ(ierr); 984 } 985 986 for (il = 0; il < Nk; il++) { 987 PetscReal val = 0.; 988 for (jl = 0; jl < dT; jl++) val += sign * wedgeMat[il * dT + jl] * workT[jl]; 989 work[il] = val; 990 } 991 if (k < 0) { 992 ierr = PetscDTAltVStar(dim, PetscAbsInt(k), -1, work, workstar);CHKERRQ(ierr); 993 #if defined(PETSC_USE_COMPLEX) 994 for (l = 0; l < Nk; l++) workS[l] = workstar[l]; 995 vals = &workS[0]; 996 #else 997 vals = &workstar[0]; 998 #endif 999 } else { 1000 #if defined(PETSC_USE_COMPLEX) 1001 for (l = 0; l < Nk; l++) workS[l] = work[l]; 1002 vals = &workS[0]; 1003 #else 1004 vals = &work[0]; 1005 #endif 1006 } 1007 for (l = 0; l < Nk; l++) { 1008 ierr = MatSetValue(prod, i, col * Nk + l, vals[l], INSERT_VALUES);CHKERRQ(ierr); 1009 } /* Nk */ 1010 } /* jT */ 1011 } /* jF */ 1012 ierr = MatRestoreRow(trace, iT, &ncolsT, &colsT, &valsT);CHKERRQ(ierr); 1013 } /* iT */ 1014 ierr = MatRestoreRow(fiber, iF, &ncolsF, &colsF, &valsF);CHKERRQ(ierr); 1015 } /* iF */ 1016 ierr = PetscFree(wedgeMat);CHKERRQ(ierr); 1017 ierr = PetscFree4(projT, projF, projTstar, projFstar);CHKERRQ(ierr); 1018 ierr = PetscFree2(workT2, workF2);CHKERRQ(ierr); 1019 ierr = PetscFree5(workT, workF, work, workstar, workS);CHKERRQ(ierr); 1020 ierr = MatAssemblyBegin(prod, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1021 ierr = MatAssemblyEnd(prod, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1022 *product = prod; 1023 PetscFunctionReturn(0); 1024 } 1025 1026 static PetscErrorCode PetscQuadraturePointsMerge(PetscQuadrature quadA, PetscQuadrature quadB, PetscQuadrature *quadJoint, PetscInt *aToJoint[], PetscInt *bToJoint[]) 1027 { 1028 PetscInt dimA, dimB; 1029 PetscInt nA, nB, nJoint, i, j, d; 1030 const PetscReal *pointsA; 1031 const PetscReal *pointsB; 1032 PetscReal *pointsJoint; 1033 PetscInt *aToJ, *bToJ; 1034 PetscQuadrature qJ; 1035 PetscErrorCode ierr; 1036 1037 PetscFunctionBegin; 1038 ierr = PetscQuadratureGetData(quadA, &dimA, NULL, &nA, &pointsA, NULL);CHKERRQ(ierr); 1039 ierr = PetscQuadratureGetData(quadB, &dimB, NULL, &nB, &pointsB, NULL);CHKERRQ(ierr); 1040 if (dimA != dimB) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Quadrature points must be in the same dimension"); 1041 nJoint = nA; 1042 ierr = PetscMalloc1(nA, &aToJ);CHKERRQ(ierr); 1043 for (i = 0; i < nA; i++) aToJ[i] = i; 1044 ierr = PetscMalloc1(nB, &bToJ);CHKERRQ(ierr); 1045 for (i = 0; i < nB; i++) { 1046 for (j = 0; j < nA; j++) { 1047 bToJ[i] = -1; 1048 for (d = 0; d < dimA; d++) if (PetscAbsReal(pointsB[i * dimA + d] - pointsA[j * dimA + d]) > PETSC_SMALL) break; 1049 if (d == dimA) { 1050 bToJ[i] = j; 1051 break; 1052 } 1053 } 1054 if (bToJ[i] == -1) { 1055 bToJ[i] = nJoint++; 1056 } 1057 } 1058 *aToJoint = aToJ; 1059 *bToJoint = bToJ; 1060 ierr = PetscMalloc1(nJoint * dimA, &pointsJoint);CHKERRQ(ierr); 1061 ierr = PetscArraycpy(pointsJoint, pointsA, nA * dimA);CHKERRQ(ierr); 1062 for (i = 0; i < nB; i++) { 1063 if (bToJ[i] >= nA) { 1064 for (d = 0; d < dimA; d++) pointsJoint[bToJ[i] * dimA + d] = pointsB[i * dimA + d]; 1065 } 1066 } 1067 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &qJ);CHKERRQ(ierr); 1068 ierr = PetscQuadratureSetData(qJ, dimA, 0, nJoint, pointsJoint, NULL);CHKERRQ(ierr); 1069 *quadJoint = qJ; 1070 PetscFunctionReturn(0); 1071 } 1072 1073 static PetscErrorCode MatricesMerge(Mat matA, Mat matB, PetscInt dim, PetscInt k, PetscInt numMerged, const PetscInt aToMerged[], const PetscInt bToMerged[], Mat *matMerged) 1074 { 1075 PetscInt m, n, mA, nA, mB, nB, Nk, i, j, l; 1076 Mat M; 1077 PetscInt *nnz; 1078 PetscInt maxnnz; 1079 PetscInt *work; 1080 PetscErrorCode ierr; 1081 1082 PetscFunctionBegin; 1083 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 1084 ierr = MatGetSize(matA, &mA, &nA);CHKERRQ(ierr); 1085 if (nA % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matA column space not a multiple of k-form size"); 1086 ierr = MatGetSize(matB, &mB, &nB);CHKERRQ(ierr); 1087 if (nB % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matB column space not a multiple of k-form size"); 1088 m = mA + mB; 1089 n = numMerged * Nk; 1090 ierr = PetscMalloc1(m, &nnz);CHKERRQ(ierr); 1091 maxnnz = 0; 1092 for (i = 0; i < mA; i++) { 1093 ierr = MatGetRow(matA, i, &(nnz[i]), NULL, NULL);CHKERRQ(ierr); 1094 if (nnz[i] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matA are not in k-form size blocks"); 1095 maxnnz = PetscMax(maxnnz, nnz[i]); 1096 } 1097 for (i = 0; i < mB; i++) { 1098 ierr = MatGetRow(matB, i, &(nnz[i+mA]), NULL, NULL);CHKERRQ(ierr); 1099 if (nnz[i+mA] % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "nonzeros in matB are not in k-form size blocks"); 1100 maxnnz = PetscMax(maxnnz, nnz[i+mA]); 1101 } 1102 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m, n, 0, nnz, &M);CHKERRQ(ierr); 1103 ierr = PetscFree(nnz);CHKERRQ(ierr); 1104 /* reasoning about which points each dof needs depends on having zeros computed at points preserved */ 1105 ierr = MatSetOption(M, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr); 1106 ierr = PetscMalloc1(maxnnz, &work);CHKERRQ(ierr); 1107 for (i = 0; i < mA; i++) { 1108 const PetscInt *cols; 1109 const PetscScalar *vals; 1110 PetscInt nCols; 1111 ierr = MatGetRow(matA, i, &nCols, &cols, &vals);CHKERRQ(ierr); 1112 for (j = 0; j < nCols / Nk; j++) { 1113 PetscInt newCol = aToMerged[cols[j * Nk] / Nk]; 1114 for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l; 1115 } 1116 ierr = MatSetValuesBlocked(M, 1, &i, nCols, work, vals, INSERT_VALUES);CHKERRQ(ierr); 1117 ierr = MatRestoreRow(matA, i, &nCols, &cols, &vals);CHKERRQ(ierr); 1118 } 1119 for (i = 0; i < mB; i++) { 1120 const PetscInt *cols; 1121 const PetscScalar *vals; 1122 1123 PetscInt row = i + mA; 1124 PetscInt nCols; 1125 ierr = MatGetRow(matB, i, &nCols, &cols, &vals);CHKERRQ(ierr); 1126 for (j = 0; j < nCols / Nk; j++) { 1127 PetscInt newCol = bToMerged[cols[j * Nk] / Nk]; 1128 for (l = 0; l < Nk; l++) work[j * Nk + l] = newCol * Nk + l; 1129 } 1130 ierr = MatSetValuesBlocked(M, 1, &row, nCols, work, vals, INSERT_VALUES);CHKERRQ(ierr); 1131 ierr = MatRestoreRow(matB, i, &nCols, &cols, &vals);CHKERRQ(ierr); 1132 } 1133 ierr = PetscFree(work);CHKERRQ(ierr); 1134 ierr = MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1135 ierr = MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1136 *matMerged = M; 1137 PetscFunctionReturn(0); 1138 } 1139 1140 static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Lagrange(PetscDualSpace sp, DM K, PetscInt f, PetscInt k, PetscInt Ncopies, PetscBool interiorOnly, PetscDualSpace *bdsp) 1141 { 1142 PetscInt Nknew, Ncnew; 1143 PetscInt dim, pointDim = -1; 1144 PetscInt depth; 1145 DM dm; 1146 PetscDualSpace_Lag *newlag; 1147 PetscErrorCode ierr; 1148 1149 PetscFunctionBegin; 1150 ierr = PetscDualSpaceGetDM(sp,&dm);CHKERRQ(ierr); 1151 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 1152 ierr = DMPlexGetDepth(dm,&depth);CHKERRQ(ierr); 1153 ierr = PetscDualSpaceDuplicate(sp,bdsp);CHKERRQ(ierr); 1154 ierr = PetscDualSpaceSetFormDegree(*bdsp,k);CHKERRQ(ierr); 1155 if (!K) { 1156 PetscBool isSimplex; 1157 1158 1159 if (depth == dim) { 1160 PetscInt coneSize; 1161 1162 pointDim = dim - 1; 1163 ierr = DMPlexGetConeSize(dm,f,&coneSize);CHKERRQ(ierr); 1164 isSimplex = (PetscBool) (coneSize == dim); 1165 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, dim-1, isSimplex, &K);CHKERRQ(ierr); 1166 } else if (depth == 1) { 1167 pointDim = 0; 1168 ierr = PetscDualSpaceCreateReferenceCell(*bdsp, 0, PETSC_TRUE, &K);CHKERRQ(ierr); 1169 } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element"); 1170 } else { 1171 ierr = PetscObjectReference((PetscObject)K);CHKERRQ(ierr); 1172 ierr = DMGetDimension(K, &pointDim);CHKERRQ(ierr); 1173 } 1174 ierr = PetscDualSpaceSetDM(*bdsp, K);CHKERRQ(ierr); 1175 ierr = DMDestroy(&K);CHKERRQ(ierr); 1176 ierr = PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew);CHKERRQ(ierr); 1177 Ncnew = Nknew * Ncopies; 1178 ierr = PetscDualSpaceSetNumComponents(*bdsp, Ncnew);CHKERRQ(ierr); 1179 newlag = (PetscDualSpace_Lag *) (*bdsp)->data; 1180 newlag->interiorOnly = interiorOnly; 1181 ierr = PetscDualSpaceSetUp(*bdsp);CHKERRQ(ierr); 1182 PetscFunctionReturn(0); 1183 } 1184 1185 static PetscErrorCode PetscDualSpaceLagrangeCreateSimplexNodeMat(Petsc1DNodeFamily nodeFamily, PetscInt dim, PetscInt sum, PetscInt Nk, PetscInt numNodeSkip, PetscQuadrature *iNodes, Mat *iMat, PetscLagNodeIndices *nodeIndices) 1186 { 1187 PetscReal *extraNodeCoords, *nodeCoords; 1188 PetscInt nNodes, nExtraNodes; 1189 PetscInt i, j, k, extraSum = sum + numNodeSkip * (1 + dim); 1190 PetscQuadrature intNodes; 1191 Mat intMat; 1192 PetscLagNodeIndices ni; 1193 PetscErrorCode ierr; 1194 1195 PetscFunctionBegin; 1196 ierr = PetscDTBinomialInt(dim + sum, dim, &nNodes);CHKERRQ(ierr); 1197 ierr = PetscDTBinomialInt(dim + extraSum, dim, &nExtraNodes);CHKERRQ(ierr); 1198 1199 ierr = PetscMalloc1(dim * nExtraNodes, &extraNodeCoords);CHKERRQ(ierr); 1200 ierr = PetscNew(&ni);CHKERRQ(ierr); 1201 ni->nodeIdxDim = dim + 1; 1202 ni->nodeVecDim = Nk; 1203 ni->nNodes = nNodes * Nk; 1204 ni->refct = 1; 1205 ierr = PetscMalloc1(nNodes * Nk * (dim + 1), &(ni->nodeIdx));CHKERRQ(ierr); 1206 ierr = PetscMalloc1(nNodes * Nk * Nk, &(ni->nodeVec));CHKERRQ(ierr); 1207 for (i = 0; i < nNodes; i++) for (j = 0; j < Nk; j++) for (k = 0; k < Nk; k++) ni->nodeVec[(i * Nk + j) * Nk + k] = (j == k) ? 1. : 0.; 1208 ierr = Petsc1DNodeFamilyComputeSimplexNodes(nodeFamily, dim, extraSum, extraNodeCoords);CHKERRQ(ierr); 1209 if (numNodeSkip) { 1210 PetscInt k; 1211 PetscInt *tup; 1212 1213 ierr = PetscMalloc1(dim * nNodes, &nodeCoords);CHKERRQ(ierr); 1214 ierr = PetscMalloc1(dim + 1, &tup);CHKERRQ(ierr); 1215 for (k = 0; k < nNodes; k++) { 1216 PetscInt j, c; 1217 PetscInt index; 1218 1219 ierr = PetscDTIndexToBary(dim + 1, sum, k, tup);CHKERRQ(ierr); 1220 for (j = 0; j < dim + 1; j++) tup[j] += numNodeSkip; 1221 for (c = 0; c < Nk; c++) { 1222 for (j = 0; j < dim + 1; j++) { 1223 ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1; 1224 } 1225 } 1226 ierr = PetscDTBaryToIndex(dim + 1, extraSum, tup, &index);CHKERRQ(ierr); 1227 for (j = 0; j < dim; j++) nodeCoords[k * dim + j] = extraNodeCoords[index * dim + j]; 1228 } 1229 ierr = PetscFree(tup);CHKERRQ(ierr); 1230 ierr = PetscFree(extraNodeCoords);CHKERRQ(ierr); 1231 } else { 1232 PetscInt k; 1233 PetscInt *tup; 1234 1235 nodeCoords = extraNodeCoords; 1236 ierr = PetscMalloc1(dim + 1, &tup);CHKERRQ(ierr); 1237 for (k = 0; k < nNodes; k++) { 1238 PetscInt j, c; 1239 1240 ierr = PetscDTIndexToBary(dim + 1, sum, k, tup);CHKERRQ(ierr); 1241 for (c = 0; c < Nk; c++) { 1242 for (j = 0; j < dim + 1; j++) { 1243 /* barycentric indices can have zeros, but we don't want to push forward zeros because it makes it harder to 1244 * determine which nodes correspond to which under symmetries, so we increase by 1 */ 1245 ni->nodeIdx[(k * Nk + c) * (dim + 1) + j] = tup[j] + 1; 1246 } 1247 } 1248 } 1249 ierr = PetscFree(tup);CHKERRQ(ierr); 1250 } 1251 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &intNodes);CHKERRQ(ierr); 1252 ierr = PetscQuadratureSetData(intNodes, dim, 0, nNodes, nodeCoords, NULL);CHKERRQ(ierr); 1253 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes * Nk, nNodes * Nk, Nk, NULL, &intMat);CHKERRQ(ierr); 1254 ierr = MatSetOption(intMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);CHKERRQ(ierr); 1255 for (j = 0; j < nNodes * Nk; j++) { 1256 PetscInt rem = j % Nk; 1257 PetscInt a, aprev = j - rem; 1258 PetscInt anext = aprev + Nk; 1259 1260 for (a = aprev; a < anext; a++) { 1261 ierr = MatSetValue(intMat, j, a, (a == j) ? 1. : 0., INSERT_VALUES);CHKERRQ(ierr); 1262 } 1263 } 1264 ierr = MatAssemblyBegin(intMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1265 ierr = MatAssemblyEnd(intMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1266 *iNodes = intNodes; 1267 *iMat = intMat; 1268 *nodeIndices = ni; 1269 PetscFunctionReturn(0); 1270 } 1271 1272 static PetscErrorCode PetscDualSpaceLagrangeCreateAllNodeIdx(PetscDualSpace sp) 1273 { 1274 DM dm; 1275 PetscInt dim, nDofs; 1276 PetscSection section; 1277 PetscInt pStart, pEnd, p; 1278 PetscInt formDegree, Nk; 1279 PetscInt nodeIdxDim, spintdim; 1280 PetscDualSpace_Lag *lag; 1281 PetscLagNodeIndices ni, verti; 1282 PetscErrorCode ierr; 1283 1284 PetscFunctionBegin; 1285 lag = (PetscDualSpace_Lag *) sp->data; 1286 verti = lag->vertIndices; 1287 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1288 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1289 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr); 1290 ierr = PetscDTBinomialInt(dim, PetscAbsInt(formDegree), &Nk);CHKERRQ(ierr); 1291 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 1292 ierr = PetscSectionGetStorageSize(section, &nDofs);CHKERRQ(ierr); 1293 ierr = PetscNew(&ni);CHKERRQ(ierr); 1294 ni->nodeIdxDim = nodeIdxDim = verti->nodeIdxDim; 1295 ni->nodeVecDim = Nk; 1296 ni->nNodes = nDofs; 1297 ni->refct = 1; 1298 ierr = PetscMalloc1(nodeIdxDim * nDofs, &(ni->nodeIdx));CHKERRQ(ierr); 1299 ierr = PetscMalloc1(Nk * nDofs, &(ni->nodeVec));CHKERRQ(ierr); 1300 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1301 ierr = PetscSectionGetDof(section, 0, &spintdim);CHKERRQ(ierr); 1302 if (spintdim) { 1303 ierr = PetscArraycpy(ni->nodeIdx, lag->intNodeIndices->nodeIdx, spintdim * nodeIdxDim);CHKERRQ(ierr); 1304 ierr = PetscArraycpy(ni->nodeVec, lag->intNodeIndices->nodeVec, spintdim * Nk);CHKERRQ(ierr); 1305 } 1306 for (p = pStart + 1; p < pEnd; p++) { 1307 PetscDualSpace psp = sp->pointSpaces[p]; 1308 PetscDualSpace_Lag *plag; 1309 PetscInt dof, off; 1310 1311 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 1312 if (!dof) continue; 1313 plag = (PetscDualSpace_Lag *) psp->data; 1314 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 1315 ierr = PetscLagNodeIndicesPushForward(dm, verti, p, plag->vertIndices, plag->intNodeIndices, 0, formDegree, &(ni->nodeIdx[off * nodeIdxDim]), &(ni->nodeVec[off * Nk]));CHKERRQ(ierr); 1316 } 1317 lag->allNodeIndices = ni; 1318 PetscFunctionReturn(0); 1319 } 1320 1321 static PetscErrorCode PetscDualSpaceCreateAllDataFromInteriorData(PetscDualSpace sp) 1322 { 1323 DM dm; 1324 PetscSection section; 1325 PetscInt pStart, pEnd, p, k, Nk, dim, Nc; 1326 PetscInt nNodes; 1327 PetscInt countNodes; 1328 Mat allMat; 1329 PetscQuadrature allNodes; 1330 PetscInt nDofs; 1331 PetscInt maxNzforms, j; 1332 PetscScalar *work; 1333 PetscReal *L, *J, *Jinv, *v0, *pv0; 1334 PetscInt *iwork; 1335 PetscReal *nodes; 1336 PetscErrorCode ierr; 1337 1338 PetscFunctionBegin; 1339 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1340 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1341 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 1342 ierr = PetscSectionGetStorageSize(section, &nDofs);CHKERRQ(ierr); 1343 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1344 ierr = PetscDualSpaceGetFormDegree(sp, &k);CHKERRQ(ierr); 1345 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 1346 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 1347 for (p = pStart, nNodes = 0, maxNzforms = 0; p < pEnd; p++) { 1348 PetscDualSpace psp; 1349 DM pdm; 1350 PetscInt pdim, pNk; 1351 PetscQuadrature intNodes; 1352 Mat intMat; 1353 1354 ierr = PetscDualSpaceGetPointSubspace(sp, p, &psp);CHKERRQ(ierr); 1355 if (!psp) continue; 1356 ierr = PetscDualSpaceGetDM(psp, &pdm);CHKERRQ(ierr); 1357 ierr = DMGetDimension(pdm, &pdim);CHKERRQ(ierr); 1358 if (pdim < PetscAbsInt(k)) continue; 1359 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);CHKERRQ(ierr); 1360 ierr = PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);CHKERRQ(ierr); 1361 if (intNodes) { 1362 PetscInt nNodesp; 1363 1364 ierr = PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, NULL, NULL);CHKERRQ(ierr); 1365 nNodes += nNodesp; 1366 } 1367 if (intMat) { 1368 PetscInt maxNzsp; 1369 PetscInt maxNzformsp; 1370 1371 ierr = MatSeqAIJGetMaxRowNonzeros(intMat, &maxNzsp);CHKERRQ(ierr); 1372 if (maxNzsp % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms"); 1373 maxNzformsp = maxNzsp / pNk; 1374 maxNzforms = PetscMax(maxNzforms, maxNzformsp); 1375 } 1376 } 1377 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nDofs, nNodes * Nc, maxNzforms * Nk, NULL, &allMat);CHKERRQ(ierr); 1378 ierr = MatSetOption(allMat,MAT_IGNORE_ZERO_ENTRIES,PETSC_FALSE);CHKERRQ(ierr); 1379 ierr = PetscMalloc7(dim, &v0, dim, &pv0, dim * dim, &J, dim * dim, &Jinv, Nk * Nk, &L, maxNzforms * Nk, &work, maxNzforms * Nk, &iwork);CHKERRQ(ierr); 1380 for (j = 0; j < dim; j++) pv0[j] = -1.; 1381 ierr = PetscMalloc1(dim * nNodes, &nodes);CHKERRQ(ierr); 1382 for (p = pStart, countNodes = 0; p < pEnd; p++) { 1383 PetscDualSpace psp; 1384 PetscQuadrature intNodes; 1385 DM pdm; 1386 PetscInt pdim, pNk; 1387 PetscInt countNodesIn = countNodes; 1388 PetscReal detJ; 1389 Mat intMat; 1390 1391 ierr = PetscDualSpaceGetPointSubspace(sp, p, &psp);CHKERRQ(ierr); 1392 if (!psp) continue; 1393 ierr = PetscDualSpaceGetDM(psp, &pdm);CHKERRQ(ierr); 1394 ierr = DMGetDimension(pdm, &pdim);CHKERRQ(ierr); 1395 if (pdim < PetscAbsInt(k)) continue; 1396 ierr = PetscDualSpaceGetInteriorData(psp, &intNodes, &intMat);CHKERRQ(ierr); 1397 if (intNodes == NULL && intMat == NULL) continue; 1398 ierr = PetscDTBinomialInt(pdim, PetscAbsInt(k), &pNk);CHKERRQ(ierr); 1399 if (p) { 1400 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, Jinv, &detJ);CHKERRQ(ierr); 1401 } else { /* identity */ 1402 PetscInt i,j; 1403 1404 for (i = 0; i < dim; i++) for (j = 0; j < dim; j++) J[i * dim + j] = Jinv[i * dim + j] = 0.; 1405 for (i = 0; i < dim; i++) J[i * dim + i] = Jinv[i * dim + i] = 1.; 1406 for (i = 0; i < dim; i++) v0[i] = -1.; 1407 } 1408 if (pdim != dim) { /* compactify Jacobian */ 1409 PetscInt i, j; 1410 1411 for (i = 0; i < dim; i++) for (j = 0; j < pdim; j++) J[i * pdim + j] = J[i * dim + j]; 1412 } 1413 ierr = PetscDTAltVPullbackMatrix(pdim, dim, J, k, L);CHKERRQ(ierr); 1414 if (intNodes) { /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L^T */ 1415 PetscInt nNodesp; 1416 const PetscReal *nodesp; 1417 PetscInt j; 1418 1419 ierr = PetscQuadratureGetData(intNodes, NULL, NULL, &nNodesp, &nodesp, NULL);CHKERRQ(ierr); 1420 for (j = 0; j < nNodesp; j++, countNodes++) { 1421 PetscInt d, e; 1422 1423 for (d = 0; d < dim; d++) { 1424 nodes[countNodes * dim + d] = v0[d]; 1425 for (e = 0; e < pdim; e++) { 1426 nodes[countNodes * dim + d] += J[d * pdim + e] * (nodesp[j * pdim + e] - pv0[e]); 1427 } 1428 } 1429 } 1430 } 1431 if (intMat) { 1432 PetscInt nrows; 1433 PetscInt off; 1434 1435 ierr = PetscSectionGetDof(section, p, &nrows);CHKERRQ(ierr); 1436 ierr = PetscSectionGetOffset(section, p, &off);CHKERRQ(ierr); 1437 for (j = 0; j < nrows; j++) { 1438 PetscInt ncols; 1439 const PetscInt *cols; 1440 const PetscScalar *vals; 1441 PetscInt l, d, e; 1442 PetscInt row = j + off; 1443 1444 ierr = MatGetRow(intMat, j, &ncols, &cols, &vals);CHKERRQ(ierr); 1445 if (ncols % pNk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms"); 1446 for (l = 0; l < ncols / pNk; l++) { 1447 PetscInt blockcol; 1448 1449 for (d = 0; d < pNk; d++) { 1450 if ((cols[l * pNk + d] % pNk) != d) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "interior matrix is not laid out as blocks of k-forms"); 1451 } 1452 blockcol = cols[l * pNk] / pNk; 1453 for (d = 0; d < Nk; d++) { 1454 iwork[l * Nk + d] = (blockcol + countNodesIn) * Nk + d; 1455 } 1456 for (d = 0; d < Nk; d++) work[l * Nk + d] = 0.; 1457 for (d = 0; d < Nk; d++) { 1458 for (e = 0; e < pNk; e++) { 1459 /* "push forward" dof by pulling back a k-form to be evaluated on the point: multiply on the right by L */ 1460 work[l * Nk + d] += vals[l * pNk + e] * L[e * pNk + d]; 1461 } 1462 } 1463 } 1464 ierr = MatSetValues(allMat, 1, &row, (ncols / pNk) * Nk, iwork, work, INSERT_VALUES);CHKERRQ(ierr); 1465 ierr = MatRestoreRow(intMat, j, &ncols, &cols, &vals);CHKERRQ(ierr); 1466 } 1467 } 1468 } 1469 ierr = MatAssemblyBegin(allMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1470 ierr = MatAssemblyEnd(allMat, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1471 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &allNodes);CHKERRQ(ierr); 1472 ierr = PetscQuadratureSetData(allNodes, dim, 0, nNodes, nodes, NULL);CHKERRQ(ierr); 1473 ierr = PetscFree7(v0, pv0, J, Jinv, L, work, iwork);CHKERRQ(ierr); 1474 ierr = MatDestroy(&(sp->allMat));CHKERRQ(ierr); 1475 sp->allMat = allMat; 1476 ierr = PetscQuadratureDestroy(&(sp->allNodes));CHKERRQ(ierr); 1477 sp->allNodes = allNodes; 1478 PetscFunctionReturn(0); 1479 } 1480 1481 static PetscErrorCode PetscDualSpaceComputeFunctionalsFromAllData(PetscDualSpace sp) 1482 { 1483 PetscQuadrature allNodes; 1484 Mat allMat; 1485 PetscInt nDofs; 1486 PetscInt dim, k, Nk, Nc, f; 1487 DM dm; 1488 PetscInt nNodes, spdim; 1489 const PetscReal *nodes = NULL; 1490 PetscSection section; 1491 PetscErrorCode ierr; 1492 1493 PetscFunctionBegin; 1494 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 1495 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1496 ierr = PetscDualSpaceGetNumComponents(sp, &Nc);CHKERRQ(ierr); 1497 ierr = PetscDualSpaceGetFormDegree(sp, &k);CHKERRQ(ierr); 1498 ierr = PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk);CHKERRQ(ierr); 1499 ierr = PetscDualSpaceGetAllData(sp, &allNodes, &allMat);CHKERRQ(ierr); 1500 nNodes = 0; 1501 if (allNodes) { 1502 ierr = PetscQuadratureGetData(allNodes, NULL, NULL, &nNodes, &nodes, NULL);CHKERRQ(ierr); 1503 } 1504 ierr = MatGetSize(allMat, &nDofs, NULL);CHKERRQ(ierr); 1505 ierr = PetscDualSpaceGetSection(sp, §ion);CHKERRQ(ierr); 1506 ierr = PetscSectionGetStorageSize(section, &spdim);CHKERRQ(ierr); 1507 if (spdim != nDofs) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "incompatible all matrix size"); 1508 ierr = PetscMalloc1(nDofs, &(sp->functional));CHKERRQ(ierr); 1509 for (f = 0; f < nDofs; f++) { 1510 PetscInt ncols, c; 1511 const PetscInt *cols; 1512 const PetscScalar *vals; 1513 PetscReal *nodesf; 1514 PetscReal *weightsf; 1515 PetscInt nNodesf; 1516 PetscInt countNodes; 1517 1518 ierr = MatGetRow(allMat, f, &ncols, &cols, &vals);CHKERRQ(ierr); 1519 if (ncols % Nk) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "all matrix is not laid out as blocks of k-forms"); 1520 for (c = 1, nNodesf = 1; c < ncols; c++) { 1521 if ((cols[c] / Nc) != (cols[c-1] / Nc)) nNodesf++; 1522 } 1523 ierr = PetscMalloc1(dim * nNodesf, &nodesf);CHKERRQ(ierr); 1524 ierr = PetscMalloc1(Nc * nNodesf, &weightsf);CHKERRQ(ierr); 1525 for (c = 0, countNodes = 0; c < ncols; c++) { 1526 if (!c || ((cols[c] / Nc) != (cols[c-1] / Nc))) { 1527 PetscInt d; 1528 1529 for (d = 0; d < Nc; d++) { 1530 weightsf[countNodes * Nc + d] = 0.; 1531 } 1532 for (d = 0; d < dim; d++) { 1533 nodesf[countNodes * dim + d] = nodes[(cols[c] / Nc) * dim + d]; 1534 } 1535 countNodes++; 1536 } 1537 weightsf[(countNodes - 1) * Nc + (cols[c] % Nc)] = PetscRealPart(vals[c]); 1538 } 1539 ierr = PetscQuadratureCreate(PETSC_COMM_SELF, &(sp->functional[f]));CHKERRQ(ierr); 1540 ierr = PetscQuadratureSetData(sp->functional[f], dim, Nc, nNodesf, nodesf, weightsf);CHKERRQ(ierr); 1541 ierr = MatRestoreRow(allMat, f, &ncols, &cols, &vals);CHKERRQ(ierr); 1542 } 1543 PetscFunctionReturn(0); 1544 } 1545 1546 /* take a matrix meant for k-forms and expand it to one for Ncopies */ 1547 static PetscErrorCode PetscDualSpaceLagrangeMatrixCreateCopies(Mat A, PetscInt Nk, PetscInt Ncopies, Mat *Abs) 1548 { 1549 PetscInt m, n, i, j, k; 1550 PetscInt maxnnz, *nnz, *iwork; 1551 Mat Ac; 1552 PetscErrorCode ierr; 1553 1554 PetscFunctionBegin; 1555 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 1556 if (n % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of columns in A %D is not a multiple of Nk %D", n, Nk); 1557 ierr = PetscMalloc1(m * Ncopies, &nnz);CHKERRQ(ierr); 1558 for (i = 0, maxnnz = 0; i < m; i++) { 1559 PetscInt innz; 1560 ierr = MatGetRow(A, i, &innz, NULL, NULL);CHKERRQ(ierr); 1561 if (innz % Nk) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "A row %D nnzs is not a multiple of Nk %D", innz, Nk); 1562 for (j = 0; j < Ncopies; j++) nnz[i * Ncopies + j] = innz; 1563 maxnnz = PetscMax(maxnnz, innz); 1564 } 1565 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, m * Ncopies, n * Ncopies, 0, nnz, &Ac);CHKERRQ(ierr); 1566 ierr = MatSetOption(Ac, MAT_IGNORE_ZERO_ENTRIES, PETSC_FALSE);CHKERRQ(ierr); 1567 ierr = PetscFree(nnz);CHKERRQ(ierr); 1568 ierr = PetscMalloc1(maxnnz, &iwork);CHKERRQ(ierr); 1569 for (i = 0; i < m; i++) { 1570 PetscInt innz; 1571 const PetscInt *cols; 1572 const PetscScalar *vals; 1573 1574 ierr = MatGetRow(A, i, &innz, &cols, &vals);CHKERRQ(ierr); 1575 for (j = 0; j < innz; j++) iwork[j] = (cols[j] / Nk) * (Nk * Ncopies) + (cols[j] % Nk); 1576 for (j = 0; j < Ncopies; j++) { 1577 PetscInt row = i * Ncopies + j; 1578 1579 ierr = MatSetValues(Ac, 1, &row, innz, iwork, vals, INSERT_VALUES);CHKERRQ(ierr); 1580 for (k = 0; k < innz; k++) iwork[k] += Nk; 1581 } 1582 ierr = MatRestoreRow(A, i, &innz, &cols, &vals);CHKERRQ(ierr); 1583 } 1584 ierr = PetscFree(iwork);CHKERRQ(ierr); 1585 ierr = MatAssemblyBegin(Ac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1586 ierr = MatAssemblyEnd(Ac, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1587 *Abs = Ac; 1588 PetscFunctionReturn(0); 1589 } 1590 1591 static PetscErrorCode DMPlexPointIsTensor_Internal_Given(DM dm, PetscInt p, PetscInt f, PetscInt f2, PetscBool *isTensor) 1592 { 1593 PetscInt coneSize, c; 1594 const PetscInt *cone; 1595 const PetscInt *fCone; 1596 const PetscInt *f2Cone; 1597 PetscInt fs[2]; 1598 PetscInt meetSize, nmeet; 1599 const PetscInt *meet; 1600 PetscErrorCode ierr; 1601 1602 PetscFunctionBegin; 1603 fs[0] = f; 1604 fs[1] = f2; 1605 ierr = DMPlexGetMeet(dm, 2, fs, &meetSize, &meet);CHKERRQ(ierr); 1606 nmeet = meetSize; 1607 ierr = DMPlexRestoreMeet(dm, 2, fs, &meetSize, &meet);CHKERRQ(ierr); 1608 if (nmeet) { 1609 *isTensor = PETSC_FALSE; 1610 PetscFunctionReturn(0); 1611 } 1612 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1613 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1614 ierr = DMPlexGetCone(dm, f, &fCone);CHKERRQ(ierr); 1615 ierr = DMPlexGetCone(dm, f2, &f2Cone);CHKERRQ(ierr); 1616 for (c = 0; c < coneSize; c++) { 1617 PetscInt e, ef; 1618 PetscInt d = -1, d2 = -1; 1619 PetscInt dcount, d2count; 1620 PetscInt t = cone[c]; 1621 PetscInt tConeSize; 1622 PetscBool tIsTensor; 1623 const PetscInt *tCone; 1624 1625 if (t == f || t == f2) continue; 1626 ierr = DMPlexGetConeSize(dm, t, &tConeSize);CHKERRQ(ierr); 1627 ierr = DMPlexGetCone(dm, t, &tCone);CHKERRQ(ierr); 1628 1629 dcount = 0; 1630 d2count = 0; 1631 for (e = 0; e < tConeSize; e++) { 1632 PetscInt q = tCone[e]; 1633 for (ef = 0; ef < coneSize - 2; ef++) { 1634 if (fCone[ef] == q) { 1635 if (dcount) { 1636 *isTensor = PETSC_FALSE; 1637 PetscFunctionReturn(0); 1638 } 1639 d = q; 1640 dcount++; 1641 } else if (f2Cone[ef] == q) { 1642 if (d2count) { 1643 *isTensor = PETSC_FALSE; 1644 PetscFunctionReturn(0); 1645 } 1646 d2 = q; 1647 d2count++; 1648 } 1649 } 1650 } 1651 ierr = DMPlexPointIsTensor_Internal_Given(dm, t, d, d2, &tIsTensor);CHKERRQ(ierr); 1652 if (!tIsTensor) { 1653 *isTensor = PETSC_FALSE; 1654 PetscFunctionReturn(0); 1655 } 1656 } 1657 *isTensor = PETSC_TRUE; 1658 PetscFunctionReturn(0); 1659 } 1660 1661 static PetscErrorCode DMPlexPointIsTensor_Internal(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB) 1662 { 1663 PetscInt coneSize, c, c2; 1664 const PetscInt *cone; 1665 PetscErrorCode ierr; 1666 1667 PetscFunctionBegin; 1668 ierr = DMPlexGetConeSize(dm, p, &coneSize);CHKERRQ(ierr); 1669 if (!coneSize) { 1670 if (isTensor) *isTensor = PETSC_FALSE; 1671 if (endA) *endA = -1; 1672 if (endB) *endB = -1; 1673 } 1674 ierr = DMPlexGetCone(dm, p, &cone);CHKERRQ(ierr); 1675 for (c = 0; c < coneSize; c++) { 1676 PetscInt f = cone[c]; 1677 PetscInt fConeSize; 1678 1679 ierr = DMPlexGetConeSize(dm, f, &fConeSize);CHKERRQ(ierr); 1680 if (fConeSize != coneSize - 2) continue; 1681 1682 for (c2 = c + 1; c2 < coneSize; c2++) { 1683 PetscInt f2 = cone[c2]; 1684 PetscBool isTensorff2; 1685 PetscInt f2ConeSize; 1686 1687 ierr = DMPlexGetConeSize(dm, f2, &f2ConeSize);CHKERRQ(ierr); 1688 if (f2ConeSize != coneSize - 2) continue; 1689 1690 ierr = DMPlexPointIsTensor_Internal_Given(dm, p, f, f2, &isTensorff2);CHKERRQ(ierr); 1691 if (isTensorff2) { 1692 if (isTensor) *isTensor = PETSC_TRUE; 1693 if (endA) *endA = f; 1694 if (endB) *endB = f2; 1695 PetscFunctionReturn(0); 1696 } 1697 } 1698 } 1699 if (isTensor) *isTensor = PETSC_FALSE; 1700 if (endA) *endA = -1; 1701 if (endB) *endB = -1; 1702 PetscFunctionReturn(0); 1703 } 1704 1705 static PetscErrorCode DMPlexPointIsTensor(DM dm, PetscInt p, PetscBool *isTensor, PetscInt *endA, PetscInt *endB) 1706 { 1707 DMPlexInterpolatedFlag interpolated; 1708 PetscErrorCode ierr; 1709 1710 PetscFunctionBegin; 1711 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1712 if (interpolated != DMPLEX_INTERPOLATED_FULL) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "Only for interpolated DMPlex's"); 1713 ierr = DMPlexPointIsTensor_Internal(dm, p, isTensor, endA, endB);CHKERRQ(ierr); 1714 PetscFunctionReturn(0); 1715 } 1716 1717 static PetscErrorCode MatPermuteByNodeIdx(Mat A, PetscLagNodeIndices ni, Mat *Aperm) 1718 { 1719 PetscInt m, n, i, j; 1720 PetscInt nodeIdxDim = ni->nodeIdxDim; 1721 PetscInt nodeVecDim = ni->nodeVecDim; 1722 PetscInt *perm; 1723 IS permIS; 1724 IS id; 1725 PetscInt *nIdxPerm; 1726 PetscReal *nVecPerm; 1727 PetscErrorCode ierr; 1728 1729 PetscFunctionBegin; 1730 ierr = PetscLagNodeIndicesGetPermutation(ni, &perm);CHKERRQ(ierr); 1731 ierr = MatGetSize(A, &m, &n);CHKERRQ(ierr); 1732 ierr = PetscMalloc1(nodeIdxDim * m, &nIdxPerm);CHKERRQ(ierr); 1733 ierr = PetscMalloc1(nodeVecDim * m, &nVecPerm);CHKERRQ(ierr); 1734 for (i = 0; i < m; i++) for (j = 0; j < nodeIdxDim; j++) nIdxPerm[i * nodeIdxDim + j] = ni->nodeIdx[perm[i] * nodeIdxDim + j]; 1735 for (i = 0; i < m; i++) for (j = 0; j < nodeVecDim; j++) nVecPerm[i * nodeVecDim + j] = ni->nodeVec[perm[i] * nodeVecDim + j]; 1736 ierr = ISCreateGeneral(PETSC_COMM_SELF, m, perm, PETSC_USE_POINTER, &permIS);CHKERRQ(ierr); 1737 ierr = ISSetPermutation(permIS);CHKERRQ(ierr); 1738 ierr = ISCreateStride(PETSC_COMM_SELF, n, 0, 1, &id);CHKERRQ(ierr); 1739 ierr = ISSetPermutation(id);CHKERRQ(ierr); 1740 ierr = MatPermute(A, permIS, id, Aperm);CHKERRQ(ierr); 1741 ierr = ISDestroy(&permIS);CHKERRQ(ierr); 1742 ierr = ISDestroy(&id);CHKERRQ(ierr); 1743 for (i = 0; i < m; i++) perm[i] = i; 1744 ierr = PetscFree(ni->nodeIdx);CHKERRQ(ierr); 1745 ierr = PetscFree(ni->nodeVec);CHKERRQ(ierr); 1746 ni->nodeIdx = nIdxPerm; 1747 ni->nodeVec = nVecPerm; 1748 PetscFunctionReturn(0); 1749 } 1750 1751 static PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp) 1752 { 1753 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 1754 DM dm = sp->dm; 1755 DM dmint = NULL; 1756 PetscInt order; 1757 PetscInt Nc = sp->Nc; 1758 MPI_Comm comm; 1759 PetscBool continuous; 1760 PetscSection section; 1761 PetscInt depth, dim, pStart, pEnd, cStart, cEnd, p, *pStratStart, *pStratEnd, d; 1762 PetscInt formDegree, Nk, Ncopies; 1763 PetscInt tensorf = -1, tensorf2 = -1; 1764 PetscBool tensorCell, tensorSpace; 1765 PetscBool uniform, trimmed; 1766 Petsc1DNodeFamily nodeFamily; 1767 PetscInt numNodeSkip; 1768 DMPlexInterpolatedFlag interpolated; 1769 PetscBool isbdm; 1770 PetscErrorCode ierr; 1771 1772 PetscFunctionBegin; 1773 /* step 1: sanitize input */ 1774 ierr = PetscObjectGetComm((PetscObject) sp, &comm);CHKERRQ(ierr); 1775 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 1776 ierr = PetscObjectTypeCompare((PetscObject)sp, "bdm", &isbdm);CHKERRQ(ierr); 1777 if (isbdm) { 1778 sp->k = -(dim-1); /* form degree of H-div */ 1779 ierr = PetscObjectChangeTypeName((PetscObject)sp, PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 1780 } 1781 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr); 1782 if (PetscAbsInt(formDegree) > dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "Form degree must be bounded by dimension"); 1783 ierr = PetscDTBinomialInt(dim,PetscAbsInt(formDegree),&Nk);CHKERRQ(ierr); 1784 if (sp->Nc <= 0 && lag->numCopies > 0) sp->Nc = Nk * lag->numCopies; 1785 Nc = sp->Nc; 1786 if (Nc % Nk) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of components is not a multiple of form degree size"); 1787 if (lag->numCopies <= 0) lag->numCopies = Nc / Nk; 1788 Ncopies = lag->numCopies; 1789 if (Nc / Nk != Ncopies) SETERRQ(comm, PETSC_ERR_ARG_INCOMP, "Number of copies * (dim choose k) != Nc"); 1790 if (!dim) sp->order = 0; 1791 order = sp->order; 1792 uniform = sp->uniform; 1793 if (!uniform) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Variable order not supported yet"); 1794 if (lag->trimmed && !formDegree) lag->trimmed = PETSC_FALSE; /* trimmed spaces are the same as full spaces for 0-forms */ 1795 if (lag->nodeType == PETSCDTNODES_DEFAULT) { 1796 lag->nodeType = PETSCDTNODES_GAUSSJACOBI; 1797 lag->nodeExponent = 0.; 1798 /* trimmed spaces don't include corner vertices, so don't use end nodes by default */ 1799 lag->endNodes = lag->trimmed ? PETSC_FALSE : PETSC_TRUE; 1800 } 1801 /* If a trimmed space and the user did choose nodes with endpoints, skip them by default */ 1802 if (lag->numNodeSkip < 0) lag->numNodeSkip = (lag->trimmed && lag->endNodes) ? 1 : 0; 1803 numNodeSkip = lag->numNodeSkip; 1804 if (lag->trimmed && !order) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zeroth order trimmed elements"); 1805 if (lag->trimmed && PetscAbsInt(formDegree) == dim) { /* convert trimmed n-forms to untrimmed of one polynomial order less */ 1806 sp->order--; 1807 order--; 1808 lag->trimmed = PETSC_FALSE; 1809 } 1810 trimmed = lag->trimmed; 1811 if (!order || PetscAbsInt(formDegree) == dim) lag->continuous = PETSC_FALSE; 1812 continuous = lag->continuous; 1813 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 1814 ierr = DMPlexGetChart(dm, &pStart, &pEnd);CHKERRQ(ierr); 1815 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 1816 if (pStart != 0 || cStart != 0) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Expect DM with chart starting at zero and cells first"); 1817 if (cEnd != 1) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Use PETSCDUALSPACEREFINED for multi-cell reference meshes"); 1818 ierr = DMPlexIsInterpolated(dm, &interpolated);CHKERRQ(ierr); 1819 if (interpolated != DMPLEX_INTERPOLATED_FULL) { 1820 ierr = DMPlexInterpolate(dm, &dmint);CHKERRQ(ierr); 1821 } else { 1822 ierr = PetscObjectReference((PetscObject)dm);CHKERRQ(ierr); 1823 dmint = dm; 1824 } 1825 tensorCell = PETSC_FALSE; 1826 if (dim > 1) { 1827 ierr = DMPlexPointIsTensor(dmint, 0, &tensorCell, &tensorf, &tensorf2);CHKERRQ(ierr); 1828 } 1829 lag->tensorCell = tensorCell; 1830 if (dim < 2 || !lag->tensorCell) lag->tensorSpace = PETSC_FALSE; 1831 tensorSpace = lag->tensorSpace; 1832 if (!lag->nodeFamily) { 1833 ierr = Petsc1DNodeFamilyCreate(lag->nodeType, lag->nodeExponent, lag->endNodes, &lag->nodeFamily);CHKERRQ(ierr); 1834 } 1835 nodeFamily = lag->nodeFamily; 1836 if (interpolated != DMPLEX_INTERPOLATED_FULL && continuous && (PetscAbsInt(formDegree) > 0 || order > 1)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Reference element won't support all boundary nodes"); 1837 1838 /* step 2: construct the boundary spaces */ 1839 ierr = PetscMalloc2(depth+1,&pStratStart,depth+1,&pStratEnd);CHKERRQ(ierr); 1840 ierr = PetscCalloc1(pEnd,&(sp->pointSpaces));CHKERRQ(ierr); 1841 for (d = 0; d <= depth; ++d) {ierr = DMPlexGetDepthStratum(dm, d, &pStratStart[d], &pStratEnd[d]);CHKERRQ(ierr);} 1842 ierr = PetscDualSpaceSectionCreate_Internal(sp, §ion);CHKERRQ(ierr); 1843 sp->pointSection = section; 1844 if (continuous && !(lag->interiorOnly)) { 1845 PetscInt h; 1846 1847 for (p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */ 1848 PetscReal v0[3]; 1849 DMPolytopeType ptype; 1850 PetscReal J[9], detJ; 1851 PetscInt q; 1852 1853 ierr = DMPlexComputeCellGeometryAffineFEM(dm, p, v0, J, NULL, &detJ);CHKERRQ(ierr); 1854 ierr = DMPlexGetCellType(dm, p, &ptype);CHKERRQ(ierr); 1855 1856 /* compare orders to previous facets: if computed, reference that dualspace */ 1857 for (q = pStratStart[depth - 1]; q < p; q++) { 1858 DMPolytopeType qtype; 1859 1860 ierr = DMPlexGetCellType(dm, q, &qtype);CHKERRQ(ierr); 1861 if (qtype == ptype) break; 1862 } 1863 if (q < p) { /* this facet has the same dual space as that one */ 1864 ierr = PetscObjectReference((PetscObject)sp->pointSpaces[q]);CHKERRQ(ierr); 1865 sp->pointSpaces[p] = sp->pointSpaces[q]; 1866 continue; 1867 } 1868 /* if not, recursively compute this dual space */ 1869 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,p,formDegree,Ncopies,PETSC_FALSE,&sp->pointSpaces[p]);CHKERRQ(ierr); 1870 } 1871 for (h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */ 1872 PetscInt hd = depth - h; 1873 PetscInt hdim = dim - h; 1874 1875 if (hdim < PetscAbsInt(formDegree)) break; 1876 for (p = pStratStart[hd]; p < pStratEnd[hd]; p++) { 1877 PetscInt suppSize, s; 1878 const PetscInt *supp; 1879 1880 ierr = DMPlexGetSupportSize(dm, p, &suppSize);CHKERRQ(ierr); 1881 ierr = DMPlexGetSupport(dm, p, &supp);CHKERRQ(ierr); 1882 for (s = 0; s < suppSize; s++) { 1883 DM qdm; 1884 PetscDualSpace qsp, psp; 1885 PetscInt c, coneSize, q; 1886 const PetscInt *cone; 1887 const PetscInt *refCone; 1888 1889 q = supp[0]; 1890 qsp = sp->pointSpaces[q]; 1891 ierr = DMPlexGetConeSize(dm, q, &coneSize);CHKERRQ(ierr); 1892 ierr = DMPlexGetCone(dm, q, &cone);CHKERRQ(ierr); 1893 for (c = 0; c < coneSize; c++) if (cone[c] == p) break; 1894 if (c == coneSize) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "cone/suppport mismatch"); 1895 ierr = PetscDualSpaceGetDM(qsp, &qdm);CHKERRQ(ierr); 1896 ierr = DMPlexGetCone(qdm, 0, &refCone);CHKERRQ(ierr); 1897 /* get the equivalent dual space from the support dual space */ 1898 ierr = PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp);CHKERRQ(ierr); 1899 if (!s) { 1900 ierr = PetscObjectReference((PetscObject)psp);CHKERRQ(ierr); 1901 sp->pointSpaces[p] = psp; 1902 } 1903 } 1904 } 1905 } 1906 for (p = 1; p < pEnd; p++) { 1907 PetscInt pspdim; 1908 if (!sp->pointSpaces[p]) continue; 1909 ierr = PetscDualSpaceGetInteriorDimension(sp->pointSpaces[p], &pspdim);CHKERRQ(ierr); 1910 ierr = PetscSectionSetDof(section, p, pspdim);CHKERRQ(ierr); 1911 } 1912 } 1913 1914 if (Ncopies > 1) { 1915 Mat intMatScalar, allMatScalar; 1916 PetscDualSpace scalarsp; 1917 PetscDualSpace_Lag *scalarlag; 1918 1919 ierr = PetscDualSpaceDuplicate(sp, &scalarsp);CHKERRQ(ierr); 1920 ierr = PetscDualSpaceSetNumComponents(scalarsp, Nk);CHKERRQ(ierr); 1921 ierr = PetscDualSpaceSetUp(scalarsp);CHKERRQ(ierr); 1922 ierr = PetscDualSpaceGetInteriorData(scalarsp, &(sp->intNodes), &intMatScalar);CHKERRQ(ierr); 1923 ierr = PetscObjectReference((PetscObject)(sp->intNodes));CHKERRQ(ierr); 1924 if (intMatScalar) {ierr = PetscDualSpaceLagrangeMatrixCreateCopies(intMatScalar, Nk, Ncopies, &(sp->intMat));CHKERRQ(ierr);} 1925 ierr = PetscDualSpaceGetAllData(scalarsp, &(sp->allNodes), &allMatScalar);CHKERRQ(ierr); 1926 ierr = PetscObjectReference((PetscObject)(sp->allNodes));CHKERRQ(ierr); 1927 ierr = PetscDualSpaceLagrangeMatrixCreateCopies(allMatScalar, Nk, Ncopies, &(sp->allMat));CHKERRQ(ierr); 1928 sp->spdim = scalarsp->spdim * Ncopies; 1929 sp->spintdim = scalarsp->spintdim * Ncopies; 1930 scalarlag = (PetscDualSpace_Lag *) scalarsp->data; 1931 ierr = PetscLagNodeIndicesReference(scalarlag->vertIndices);CHKERRQ(ierr); 1932 lag->vertIndices = scalarlag->vertIndices; 1933 ierr = PetscLagNodeIndicesReference(scalarlag->intNodeIndices);CHKERRQ(ierr); 1934 lag->intNodeIndices = scalarlag->intNodeIndices; 1935 ierr = PetscLagNodeIndicesReference(scalarlag->allNodeIndices);CHKERRQ(ierr); 1936 lag->allNodeIndices = scalarlag->allNodeIndices; 1937 ierr = PetscDualSpaceDestroy(&scalarsp);CHKERRQ(ierr); 1938 ierr = PetscSectionSetDof(section, 0, sp->spintdim);CHKERRQ(ierr); 1939 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 1940 ierr = PetscDualSpaceComputeFunctionalsFromAllData(sp);CHKERRQ(ierr); 1941 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr); 1942 ierr = DMDestroy(&dmint);CHKERRQ(ierr); 1943 PetscFunctionReturn(0); 1944 } 1945 1946 if (trimmed && !continuous) { 1947 /* the dofs of a trimmed space don't have a nice tensor/lattice structure: 1948 * just construct the continuous dual space and copy all of the data over, 1949 * allocating it all to the cell instead of splitting it up between the boundaries */ 1950 PetscDualSpace spcont; 1951 PetscInt spdim, f; 1952 PetscQuadrature allNodes; 1953 PetscDualSpace_Lag *lagc; 1954 Mat allMat; 1955 1956 ierr = PetscDualSpaceDuplicate(sp, &spcont);CHKERRQ(ierr); 1957 ierr = PetscDualSpaceLagrangeSetContinuity(spcont, PETSC_TRUE);CHKERRQ(ierr); 1958 ierr = PetscDualSpaceSetUp(spcont);CHKERRQ(ierr); 1959 ierr = PetscDualSpaceGetDimension(spcont, &spdim);CHKERRQ(ierr); 1960 sp->spdim = sp->spintdim = spdim; 1961 ierr = PetscSectionSetDof(section, 0, spdim);CHKERRQ(ierr); 1962 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 1963 ierr = PetscMalloc1(spdim, &(sp->functional));CHKERRQ(ierr); 1964 for (f = 0; f < spdim; f++) { 1965 PetscQuadrature fn; 1966 1967 ierr = PetscDualSpaceGetFunctional(spcont, f, &fn);CHKERRQ(ierr); 1968 ierr = PetscObjectReference((PetscObject)fn);CHKERRQ(ierr); 1969 sp->functional[f] = fn; 1970 } 1971 ierr = PetscDualSpaceGetAllData(spcont, &allNodes, &allMat);CHKERRQ(ierr); 1972 ierr = PetscObjectReference((PetscObject) allNodes);CHKERRQ(ierr); 1973 ierr = PetscObjectReference((PetscObject) allNodes);CHKERRQ(ierr); 1974 sp->allNodes = sp->intNodes = allNodes; 1975 ierr = PetscObjectReference((PetscObject) allMat);CHKERRQ(ierr); 1976 ierr = PetscObjectReference((PetscObject) allMat);CHKERRQ(ierr); 1977 sp->allMat = sp->intMat = allMat; 1978 /* TODO: copy over symmetries */ 1979 lagc = (PetscDualSpace_Lag *) spcont->data; 1980 ierr = PetscLagNodeIndicesReference(lagc->vertIndices);CHKERRQ(ierr); 1981 lag->vertIndices = lagc->vertIndices; 1982 ierr = PetscLagNodeIndicesReference(lagc->allNodeIndices);CHKERRQ(ierr); 1983 ierr = PetscLagNodeIndicesReference(lagc->allNodeIndices);CHKERRQ(ierr); 1984 lag->intNodeIndices = lagc->allNodeIndices; 1985 lag->allNodeIndices = lagc->allNodeIndices; 1986 ierr = PetscDualSpaceDestroy(&spcont);CHKERRQ(ierr); 1987 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr); 1988 ierr = DMDestroy(&dmint);CHKERRQ(ierr); 1989 PetscFunctionReturn(0); 1990 } 1991 1992 /* step 3: construct intNodes, and intMat, and combine it with boundray data to make allNodes and allMat */ 1993 if (!tensorSpace) { 1994 if (!tensorCell) {ierr = PetscLagNodeIndicesCreateSimplexVertices(dm, &(lag->vertIndices));CHKERRQ(ierr);} 1995 1996 if (trimmed) { 1997 if (order + PetscAbsInt(formDegree) > dim) { 1998 PetscInt sum = order + PetscAbsInt(formDegree) - dim - 1; 1999 PetscInt nDofs; 2000 2001 ierr = PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));CHKERRQ(ierr); 2002 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr); 2003 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr); 2004 } 2005 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2006 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr); 2007 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr); 2008 } else { 2009 if (!continuous) { 2010 PetscInt sum = order; 2011 PetscInt nDofs; 2012 2013 ierr = PetscDualSpaceLagrangeCreateSimplexNodeMat(nodeFamily, dim, sum, Nk, numNodeSkip, &sp->intNodes, &sp->intMat, &(lag->intNodeIndices));CHKERRQ(ierr); 2014 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr); 2015 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr); 2016 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2017 ierr = PetscObjectReference((PetscObject)(sp->intNodes));CHKERRQ(ierr); 2018 sp->allNodes = sp->intNodes; 2019 ierr = PetscObjectReference((PetscObject)(sp->intMat));CHKERRQ(ierr); 2020 sp->allMat = sp->intMat; 2021 ierr = PetscLagNodeIndicesReference(lag->intNodeIndices);CHKERRQ(ierr); 2022 lag->allNodeIndices = lag->intNodeIndices; 2023 } else { 2024 if (order + PetscAbsInt(formDegree) > dim) { 2025 PetscDualSpace trimmedsp; 2026 PetscDualSpace_Lag *trimmedlag; 2027 PetscQuadrature intNodes; 2028 PetscInt trFormDegree = formDegree >= 0 ? formDegree - dim : dim - PetscAbsInt(formDegree); 2029 PetscInt nDofs; 2030 Mat intMat; 2031 2032 ierr = PetscDualSpaceDuplicate(sp, &trimmedsp);CHKERRQ(ierr); 2033 ierr = PetscDualSpaceLagrangeSetTrimmed(trimmedsp, PETSC_TRUE);CHKERRQ(ierr); 2034 ierr = PetscDualSpaceSetOrder(trimmedsp, order + PetscAbsInt(formDegree) - dim);CHKERRQ(ierr); 2035 ierr = PetscDualSpaceSetFormDegree(trimmedsp, trFormDegree);CHKERRQ(ierr); 2036 trimmedlag = (PetscDualSpace_Lag *) trimmedsp->data; 2037 trimmedlag->numNodeSkip = numNodeSkip + 1; 2038 ierr = PetscDualSpaceSetUp(trimmedsp);CHKERRQ(ierr); 2039 ierr = PetscDualSpaceGetAllData(trimmedsp, &intNodes, &intMat);CHKERRQ(ierr); 2040 ierr = PetscObjectReference((PetscObject)intNodes);CHKERRQ(ierr); 2041 sp->intNodes = intNodes; 2042 ierr = PetscObjectReference((PetscObject)intMat);CHKERRQ(ierr); 2043 sp->intMat = intMat; 2044 ierr = MatGetSize(sp->intMat, &nDofs, NULL);CHKERRQ(ierr); 2045 ierr = PetscLagNodeIndicesReference(trimmedlag->allNodeIndices);CHKERRQ(ierr); 2046 lag->intNodeIndices = trimmedlag->allNodeIndices; 2047 ierr = PetscDualSpaceDestroy(&trimmedsp);CHKERRQ(ierr); 2048 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr); 2049 } 2050 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2051 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr); 2052 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr); 2053 } 2054 } 2055 } else { 2056 /* assume the tensor element has the first facet being the cross-section, having its normal 2057 * pointing in the last coordinate direction */ 2058 PetscQuadrature intNodesTrace = NULL; 2059 PetscQuadrature intNodesFiber = NULL; 2060 PetscQuadrature intNodes = NULL; 2061 PetscLagNodeIndices intNodeIndices = NULL; 2062 Mat intMat = NULL; 2063 2064 if (PetscAbsInt(formDegree) < dim) { /* get the trace k-forms on the first facet, and the 0-forms on the edge */ 2065 PetscDualSpace trace, fiber; 2066 PetscDualSpace_Lag *tracel, *fiberl; 2067 Mat intMatTrace, intMatFiber; 2068 2069 if (sp->pointSpaces[tensorf]) { 2070 ierr = PetscObjectReference((PetscObject)(sp->pointSpaces[tensorf]));CHKERRQ(ierr); 2071 trace = sp->pointSpaces[tensorf]; 2072 } else { 2073 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,formDegree,Ncopies,PETSC_TRUE,&trace);CHKERRQ(ierr); 2074 } 2075 ierr = PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,0,1,PETSC_TRUE,&fiber);CHKERRQ(ierr); 2076 tracel = (PetscDualSpace_Lag *) trace->data; 2077 fiberl = (PetscDualSpace_Lag *) fiber->data; 2078 ierr = PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));CHKERRQ(ierr); 2079 ierr = PetscDualSpaceGetInteriorData(trace, &intNodesTrace, &intMatTrace);CHKERRQ(ierr); 2080 ierr = PetscDualSpaceGetInteriorData(fiber, &intNodesFiber, &intMatFiber);CHKERRQ(ierr); 2081 if (intNodesTrace && intNodesFiber) { 2082 ierr = PetscQuadratureCreateTensor(intNodesTrace, intNodesFiber, &intNodes);CHKERRQ(ierr); 2083 ierr = MatTensorAltV(intMatTrace, intMatFiber, dim-1, formDegree, 1, 0, &intMat);CHKERRQ(ierr); 2084 ierr = PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, formDegree, fiberl->intNodeIndices, 1, 0, &intNodeIndices);CHKERRQ(ierr); 2085 } 2086 ierr = PetscObjectReference((PetscObject) intNodesTrace);CHKERRQ(ierr); 2087 ierr = PetscObjectReference((PetscObject) intNodesFiber);CHKERRQ(ierr); 2088 ierr = PetscDualSpaceDestroy(&fiber);CHKERRQ(ierr); 2089 ierr = PetscDualSpaceDestroy(&trace);CHKERRQ(ierr); 2090 } 2091 if (PetscAbsInt(formDegree) > 0) { /* get the trace (k-1)-forms on the first facet, and the 1-forms on the edge */ 2092 PetscDualSpace trace, fiber; 2093 PetscDualSpace_Lag *tracel, *fiberl; 2094 PetscQuadrature intNodesTrace2, intNodesFiber2, intNodes2; 2095 PetscLagNodeIndices intNodeIndices2; 2096 Mat intMatTrace, intMatFiber, intMat2; 2097 PetscInt traceDegree = formDegree > 0 ? formDegree - 1 : formDegree + 1; 2098 PetscInt fiberDegree = formDegree > 0 ? 1 : -1; 2099 2100 ierr = PetscDualSpaceCreateFacetSubspace_Lagrange(sp,NULL,tensorf,traceDegree,Ncopies,PETSC_TRUE,&trace);CHKERRQ(ierr); 2101 ierr = PetscDualSpaceCreateEdgeSubspace_Lagrange(sp,order,fiberDegree,1,PETSC_TRUE,&fiber);CHKERRQ(ierr); 2102 tracel = (PetscDualSpace_Lag *) trace->data; 2103 fiberl = (PetscDualSpace_Lag *) fiber->data; 2104 if (!lag->vertIndices) { 2105 ierr = PetscLagNodeIndicesCreateTensorVertices(dm, tracel->vertIndices, &(lag->vertIndices));CHKERRQ(ierr); 2106 } 2107 ierr = PetscDualSpaceGetInteriorData(trace, &intNodesTrace2, &intMatTrace);CHKERRQ(ierr); 2108 ierr = PetscDualSpaceGetInteriorData(fiber, &intNodesFiber2, &intMatFiber);CHKERRQ(ierr); 2109 if (intNodesTrace2 && intNodesFiber2) { 2110 ierr = PetscQuadratureCreateTensor(intNodesTrace2, intNodesFiber2, &intNodes2);CHKERRQ(ierr); 2111 ierr = MatTensorAltV(intMatTrace, intMatFiber, dim-1, traceDegree, 1, fiberDegree, &intMat2);CHKERRQ(ierr); 2112 ierr = PetscLagNodeIndicesTensor(tracel->intNodeIndices, dim - 1, traceDegree, fiberl->intNodeIndices, 1, fiberDegree, &intNodeIndices2);CHKERRQ(ierr); 2113 if (!intMat) { 2114 intMat = intMat2; 2115 intNodes = intNodes2; 2116 intNodeIndices = intNodeIndices2; 2117 } else { 2118 /* merge the two matrices and the two sets of points */ 2119 PetscInt nM; 2120 PetscInt nDof, nDof2; 2121 PetscInt *toMerged = NULL, *toMerged2 = NULL; 2122 PetscQuadrature merged = NULL; 2123 PetscLagNodeIndices intNodeIndicesMerged = NULL; 2124 Mat matMerged = NULL; 2125 2126 ierr = MatGetSize(intMat, &nDof, 0);CHKERRQ(ierr); 2127 ierr = MatGetSize(intMat2, &nDof2, 0);CHKERRQ(ierr); 2128 ierr = PetscQuadraturePointsMerge(intNodes, intNodes2, &merged, &toMerged, &toMerged2);CHKERRQ(ierr); 2129 ierr = PetscQuadratureGetData(merged, NULL, NULL, &nM, NULL, NULL);CHKERRQ(ierr); 2130 ierr = MatricesMerge(intMat, intMat2, dim, formDegree, nM, toMerged, toMerged2, &matMerged);CHKERRQ(ierr); 2131 ierr = PetscLagNodeIndicesMerge(intNodeIndices, intNodeIndices2, &intNodeIndicesMerged);CHKERRQ(ierr); 2132 ierr = PetscFree(toMerged);CHKERRQ(ierr); 2133 ierr = PetscFree(toMerged2);CHKERRQ(ierr); 2134 ierr = MatDestroy(&intMat);CHKERRQ(ierr); 2135 ierr = MatDestroy(&intMat2);CHKERRQ(ierr); 2136 ierr = PetscQuadratureDestroy(&intNodes);CHKERRQ(ierr); 2137 ierr = PetscQuadratureDestroy(&intNodes2);CHKERRQ(ierr); 2138 ierr = PetscLagNodeIndicesDestroy(&intNodeIndices);CHKERRQ(ierr); 2139 ierr = PetscLagNodeIndicesDestroy(&intNodeIndices2);CHKERRQ(ierr); 2140 intNodes = merged; 2141 intMat = matMerged; 2142 intNodeIndices = intNodeIndicesMerged; 2143 if (!trimmed) { 2144 Mat intMatPerm; 2145 2146 ierr = MatPermuteByNodeIdx(intMat, intNodeIndices, &intMatPerm);CHKERRQ(ierr); 2147 ierr = MatDestroy(&intMat);CHKERRQ(ierr); 2148 intMat = intMatPerm; 2149 } 2150 } 2151 } 2152 ierr = PetscDualSpaceDestroy(&fiber);CHKERRQ(ierr); 2153 ierr = PetscDualSpaceDestroy(&trace);CHKERRQ(ierr); 2154 } 2155 ierr = PetscQuadratureDestroy(&intNodesTrace);CHKERRQ(ierr); 2156 ierr = PetscQuadratureDestroy(&intNodesFiber);CHKERRQ(ierr); 2157 sp->intNodes = intNodes; 2158 sp->intMat = intMat; 2159 lag->intNodeIndices = intNodeIndices; 2160 { 2161 PetscInt nDofs = 0; 2162 2163 if (intMat) { 2164 ierr = MatGetSize(intMat, &nDofs, NULL);CHKERRQ(ierr); 2165 } 2166 ierr = PetscSectionSetDof(section, 0, nDofs);CHKERRQ(ierr); 2167 } 2168 ierr = PetscDualSpaceSectionSetUp_Internal(sp, section);CHKERRQ(ierr); 2169 if (continuous) { 2170 ierr = PetscDualSpaceCreateAllDataFromInteriorData(sp);CHKERRQ(ierr); 2171 ierr = PetscDualSpaceLagrangeCreateAllNodeIdx(sp);CHKERRQ(ierr); 2172 } else { 2173 ierr = PetscObjectReference((PetscObject) intNodes);CHKERRQ(ierr); 2174 sp->allNodes = intNodes; 2175 ierr = PetscObjectReference((PetscObject) intMat);CHKERRQ(ierr); 2176 sp->allMat = intMat; 2177 ierr = PetscLagNodeIndicesReference(intNodeIndices);CHKERRQ(ierr); 2178 lag->allNodeIndices = intNodeIndices; 2179 } 2180 } 2181 ierr = PetscSectionGetStorageSize(section, &sp->spdim);CHKERRQ(ierr); 2182 ierr = PetscSectionGetConstrainedStorageSize(section, &sp->spintdim);CHKERRQ(ierr); 2183 ierr = PetscDualSpaceComputeFunctionalsFromAllData(sp);CHKERRQ(ierr); 2184 ierr = PetscFree2(pStratStart, pStratEnd);CHKERRQ(ierr); 2185 ierr = DMDestroy(&dmint);CHKERRQ(ierr); 2186 PetscFunctionReturn(0); 2187 } 2188 2189 PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat) 2190 { 2191 PetscDualSpace_Lag *lag; 2192 DM dm; 2193 PetscLagNodeIndices vertIndices, intNodeIndices; 2194 PetscLagNodeIndices ni; 2195 PetscInt nodeIdxDim, nodeVecDim, nNodes; 2196 PetscInt formDegree; 2197 PetscInt *perm, *permOrnt; 2198 PetscInt *nnz; 2199 PetscInt n; 2200 PetscInt maxGroupSize; 2201 PetscScalar *V, *W, *work; 2202 Mat A; 2203 PetscErrorCode ierr; 2204 2205 PetscFunctionBegin; 2206 if (!sp->spintdim) { 2207 *symMat = NULL; 2208 PetscFunctionReturn(0); 2209 } 2210 lag = (PetscDualSpace_Lag *) sp->data; 2211 vertIndices = lag->vertIndices; 2212 intNodeIndices = lag->intNodeIndices; 2213 ierr = PetscDualSpaceGetDM(sp, &dm);CHKERRQ(ierr); 2214 ierr = PetscDualSpaceGetFormDegree(sp, &formDegree);CHKERRQ(ierr); 2215 ierr = PetscNew(&ni);CHKERRQ(ierr); 2216 ni->refct = 1; 2217 ni->nodeIdxDim = nodeIdxDim = intNodeIndices->nodeIdxDim; 2218 ni->nodeVecDim = nodeVecDim = intNodeIndices->nodeVecDim; 2219 ni->nNodes = nNodes = intNodeIndices->nNodes; 2220 ierr = PetscMalloc1(nNodes * nodeIdxDim, &(ni->nodeIdx));CHKERRQ(ierr); 2221 ierr = PetscMalloc1(nNodes * nodeVecDim, &(ni->nodeVec));CHKERRQ(ierr); 2222 ierr = PetscLagNodeIndicesPushForward(dm, vertIndices, 0, vertIndices, intNodeIndices, ornt, formDegree, ni->nodeIdx, ni->nodeVec);CHKERRQ(ierr); 2223 ierr = PetscLagNodeIndicesGetPermutation(intNodeIndices, &perm);CHKERRQ(ierr); 2224 ierr = PetscLagNodeIndicesGetPermutation(ni, &permOrnt);CHKERRQ(ierr); 2225 ierr = PetscMalloc1(nNodes, &nnz);CHKERRQ(ierr); 2226 for (n = 0, maxGroupSize = 0; n < nNodes;) { /* incremented in the loop */ 2227 PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]); 2228 PetscInt m, nEnd; 2229 PetscInt groupSize; 2230 for (nEnd = n + 1; nEnd < nNodes; nEnd++) { 2231 PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]); 2232 PetscInt d; 2233 2234 /* compare the oriented permutation indices */ 2235 for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break; 2236 if (d < nodeIdxDim) break; 2237 } 2238 #if defined(PETSC_USE_DEBUG) 2239 { 2240 PetscInt m; 2241 PetscInt *nind = &(intNodeIndices->nodeIdx[perm[n] * nodeIdxDim]); 2242 2243 for (m = n + 1; m < nEnd; m++) { 2244 PetscInt *mind = &(intNodeIndices->nodeIdx[perm[m] * nodeIdxDim]); 2245 PetscInt d; 2246 2247 /* compare the oriented permutation indices */ 2248 for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break; 2249 if (d < nodeIdxDim) break; 2250 } 2251 if (m < nEnd) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs with same index after symmetry not same block size"); 2252 } 2253 #endif 2254 groupSize = nEnd - n; 2255 for (m = n; m < nEnd; m++) nnz[permOrnt[m]] = groupSize; 2256 2257 maxGroupSize = PetscMax(maxGroupSize, nEnd - n); 2258 /* permOrnt[[n, nEnd)] is a group of dofs that, under the symmetry are at the same location */ 2259 n = nEnd; 2260 } 2261 if (maxGroupSize > nodeVecDim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Dofs are not in blocks that can be solved"); 2262 ierr = MatCreateSeqAIJ(PETSC_COMM_SELF, nNodes, nNodes, 0, nnz, &A);CHKERRQ(ierr); 2263 ierr = PetscFree(nnz);CHKERRQ(ierr); 2264 ierr = PetscMalloc3(maxGroupSize * nodeVecDim, &V, maxGroupSize * nodeVecDim, &W, nodeVecDim * 2, &work);CHKERRQ(ierr); 2265 for (n = 0; n < nNodes;) { /* incremented in the loop */ 2266 PetscInt *nind = &(ni->nodeIdx[permOrnt[n] * nodeIdxDim]); 2267 PetscInt nEnd; 2268 PetscInt m; 2269 PetscInt groupSize; 2270 for (nEnd = n + 1; nEnd < nNodes; nEnd++) { 2271 PetscInt *mind = &(ni->nodeIdx[permOrnt[nEnd] * nodeIdxDim]); 2272 PetscInt d; 2273 2274 /* compare the oriented permutation indices */ 2275 for (d = 0; d < nodeIdxDim; d++) if (mind[d] != nind[d]) break; 2276 if (d < nodeIdxDim) break; 2277 } 2278 groupSize = nEnd - n; 2279 /* get all of the vectors from the original */ 2280 for (m = n; m < nEnd; m++) { 2281 PetscInt d; 2282 2283 for (d = 0; d < nodeVecDim; d++) { 2284 V[(m - n) * nodeVecDim + d] = intNodeIndices->nodeVec[perm[m] * nodeVecDim + d]; 2285 W[(m - n) * nodeVecDim + d] = ni->nodeVec[permOrnt[m] * nodeVecDim + d]; 2286 } 2287 } 2288 /* now we have to solve for W in terms of V */ 2289 { 2290 char transpose = 'N'; 2291 PetscBLASInt bm = nodeVecDim; 2292 PetscBLASInt bn = groupSize; 2293 PetscBLASInt bnrhs = groupSize; 2294 PetscBLASInt blda = bm; 2295 PetscBLASInt bldb = bm; 2296 PetscBLASInt blwork = 2 * nodeVecDim; 2297 PetscBLASInt info; 2298 2299 PetscStackCallBLAS("LAPACKgels",LAPACKgels_(&transpose,&bm,&bn,&bnrhs,V,&blda,W,&bldb,work,&blwork, &info)); 2300 if (info != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Bad argument to GELS"); 2301 /* repack */ 2302 { 2303 PetscInt i, j; 2304 2305 for (i = 0; i < groupSize; i++) { 2306 for (j = 0; j < groupSize; j++) { 2307 V[i * groupSize + j] = W[i * nodeVecDim + j]; 2308 } 2309 } 2310 } 2311 } 2312 ierr = MatSetValues(A, groupSize, &permOrnt[n], groupSize, &perm[n], V, INSERT_VALUES);CHKERRQ(ierr); 2313 /* permOrnt[[n, nEnd)] is a group of dofs that, under the symmetry are at the same location */ 2314 n = nEnd; 2315 } 2316 ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2317 ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2318 *symMat = A; 2319 ierr = PetscFree3(V,W,work);CHKERRQ(ierr); 2320 ierr = PetscLagNodeIndicesDestroy(&ni);CHKERRQ(ierr); 2321 PetscFunctionReturn(0); 2322 } 2323 2324 #define BaryIndex(perEdge,a,b,c) (((b)*(2*perEdge+1-(b)))/2)+(c) 2325 2326 #define CartIndex(perEdge,a,b) (perEdge*(a)+b) 2327 2328 static PetscErrorCode PetscDualSpaceGetSymmetries_Lagrange(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips) 2329 { 2330 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 2331 PetscInt dim, order, Nc; 2332 PetscErrorCode ierr; 2333 2334 PetscFunctionBegin; 2335 ierr = PetscDualSpaceGetOrder(sp,&order);CHKERRQ(ierr); 2336 ierr = PetscDualSpaceGetNumComponents(sp,&Nc);CHKERRQ(ierr); 2337 ierr = DMGetDimension(sp->dm,&dim);CHKERRQ(ierr); 2338 if (!lag->symComputed) { /* store symmetries */ 2339 PetscInt pStart, pEnd, p; 2340 PetscInt numPoints; 2341 PetscInt numFaces; 2342 PetscInt spintdim; 2343 PetscInt ***symperms; 2344 PetscScalar ***symflips; 2345 2346 ierr = DMPlexGetChart(sp->dm, &pStart, &pEnd);CHKERRQ(ierr); 2347 numPoints = pEnd - pStart; 2348 ierr = DMPlexGetConeSize(sp->dm, 0, &numFaces);CHKERRQ(ierr); 2349 ierr = PetscCalloc1(numPoints,&symperms);CHKERRQ(ierr); 2350 ierr = PetscCalloc1(numPoints,&symflips);CHKERRQ(ierr); 2351 spintdim = sp->spintdim; 2352 /* The nodal symmetry behavior is not present when tensorSpace != tensorCell: someone might want this for the "S" 2353 * family of FEEC spaces. Most used in particular are discontinuous polynomial L2 spaces in tensor cells, where 2354 * the symmetries are not necessary for FE assembly. So for now we assume this is the case and don't return 2355 * symmetries if tensorSpace != tensorCell */ 2356 if (spintdim && 0 < dim && dim < 3 && (lag->tensorSpace == lag->tensorCell)) { /* compute self symmetries */ 2357 PetscInt **cellSymperms; 2358 PetscScalar **cellSymflips; 2359 PetscInt ornt; 2360 PetscInt nCopies = Nc / lag->intNodeIndices->nodeVecDim; 2361 PetscInt nNodes = lag->intNodeIndices->nNodes; 2362 2363 lag->numSelfSym = 2 * numFaces; 2364 lag->selfSymOff = numFaces; 2365 ierr = PetscCalloc1(2*numFaces,&cellSymperms);CHKERRQ(ierr); 2366 ierr = PetscCalloc1(2*numFaces,&cellSymflips);CHKERRQ(ierr); 2367 /* we want to be able to index symmetries directly with the orientations, which range from [-numFaces,numFaces) */ 2368 symperms[0] = &cellSymperms[numFaces]; 2369 symflips[0] = &cellSymflips[numFaces]; 2370 if (lag->intNodeIndices->nodeVecDim * nCopies != Nc) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs"); 2371 if (nNodes * nCopies != spintdim) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Node indices incompatible with dofs"); 2372 for (ornt = -numFaces; ornt < numFaces; ornt++) { /* for every symmetry, compute the symmetry matrix, and extract rows to see if it fits in the perm + flip framework */ 2373 Mat symMat; 2374 PetscInt *perm; 2375 PetscScalar *flips; 2376 PetscInt i; 2377 2378 if (!ornt) continue; 2379 ierr = PetscMalloc1(spintdim, &perm);CHKERRQ(ierr); 2380 ierr = PetscCalloc1(spintdim, &flips);CHKERRQ(ierr); 2381 for (i = 0; i < spintdim; i++) perm[i] = -1; 2382 ierr = PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(sp, ornt, &symMat);CHKERRQ(ierr); 2383 for (i = 0; i < nNodes; i++) { 2384 PetscInt ncols; 2385 PetscInt j, k; 2386 const PetscInt *cols; 2387 const PetscScalar *vals; 2388 PetscBool nz_seen = PETSC_FALSE; 2389 2390 ierr = MatGetRow(symMat, i, &ncols, &cols, &vals);CHKERRQ(ierr); 2391 for (j = 0; j < ncols; j++) { 2392 if (PetscAbsScalar(vals[j]) > PETSC_SMALL) { 2393 if (nz_seen) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips"); 2394 nz_seen = PETSC_TRUE; 2395 if (PetscAbsScalar(PetscAbsScalar(vals[j]) - 1.) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips"); 2396 if (PetscAbsReal(PetscImaginaryPart(vals[j])) > PETSC_SMALL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips"); 2397 if (perm[cols[j] * nCopies] >= 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "This dual space has symmetries that can't be described as a permutation + sign flips"); 2398 for (k = 0; k < nCopies; k++) { 2399 perm[cols[j] * nCopies + k] = i * nCopies + k; 2400 } 2401 if (PetscRealPart(vals[j]) < 0.) { 2402 for (k = 0; k < nCopies; k++) { 2403 flips[i * nCopies + k] = -1.; 2404 } 2405 } else { 2406 for (k = 0; k < nCopies; k++) { 2407 flips[i * nCopies + k] = 1.; 2408 } 2409 } 2410 } 2411 } 2412 ierr = MatRestoreRow(symMat, i, &ncols, &cols, &vals);CHKERRQ(ierr); 2413 } 2414 ierr = MatDestroy(&symMat);CHKERRQ(ierr); 2415 /* if there were no sign flips, keep NULL */ 2416 for (i = 0; i < spintdim; i++) if (flips[i] != 1.) break; 2417 if (i == spintdim) { 2418 ierr = PetscFree(flips);CHKERRQ(ierr); 2419 flips = NULL; 2420 } 2421 /* if the permutation is identity, keep NULL */ 2422 for (i = 0; i < spintdim; i++) if (perm[i] != i) break; 2423 if (i == spintdim) { 2424 ierr = PetscFree(perm);CHKERRQ(ierr); 2425 perm = NULL; 2426 } 2427 symperms[0][ornt] = perm; 2428 symflips[0][ornt] = flips; 2429 } 2430 /* if no orientations produced non-identity permutations, keep NULL */ 2431 for (ornt = -numFaces; ornt < numFaces; ornt++) if (symperms[0][ornt]) break; 2432 if (ornt == numFaces) { 2433 ierr = PetscFree(cellSymperms);CHKERRQ(ierr); 2434 symperms[0] = NULL; 2435 } 2436 /* if no orientations produced sign flips, keep NULL */ 2437 for (ornt = -numFaces; ornt < numFaces; ornt++) if (symflips[0][ornt]) break; 2438 if (ornt == numFaces) { 2439 ierr = PetscFree(cellSymflips);CHKERRQ(ierr); 2440 symflips[0] = NULL; 2441 } 2442 } 2443 { 2444 PetscInt closureSize = 0; 2445 PetscInt *closure = NULL; 2446 PetscInt r; 2447 2448 ierr = DMPlexGetTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2449 for (r = 0; r < closureSize; r++) { 2450 PetscDualSpace psp; 2451 PetscInt point = closure[2 * r]; 2452 PetscInt pspintdim; 2453 const PetscInt ***psymperms = NULL; 2454 const PetscScalar ***psymflips = NULL; 2455 2456 if (!point) continue; 2457 ierr = PetscDualSpaceGetPointSubspace(sp, point, &psp);CHKERRQ(ierr); 2458 if (!psp) continue; 2459 ierr = PetscDualSpaceGetInteriorDimension(psp, &pspintdim);CHKERRQ(ierr); 2460 if (!pspintdim) continue; 2461 ierr = PetscDualSpaceGetSymmetries(psp,&psymperms,&psymflips);CHKERRQ(ierr); 2462 symperms[r] = (PetscInt **) (psymperms ? psymperms[0] : NULL); 2463 symflips[r] = (PetscScalar **) (psymflips ? psymflips[0] : NULL); 2464 } 2465 ierr = DMPlexRestoreTransitiveClosure(sp->dm,0,PETSC_TRUE,&closureSize,&closure);CHKERRQ(ierr); 2466 } 2467 for (p = 0; p < pEnd; p++) if (symperms[p]) break; 2468 if (p == pEnd) { 2469 ierr = PetscFree(symperms);CHKERRQ(ierr); 2470 symperms = NULL; 2471 } 2472 for (p = 0; p < pEnd; p++) if (symflips[p]) break; 2473 if (p == pEnd) { 2474 ierr = PetscFree(symflips);CHKERRQ(ierr); 2475 symflips = NULL; 2476 } 2477 lag->symperms = symperms; 2478 lag->symflips = symflips; 2479 lag->symComputed = PETSC_TRUE; 2480 } 2481 if (perms) *perms = (const PetscInt ***) lag->symperms; 2482 if (flips) *flips = (const PetscScalar ***) lag->symflips; 2483 PetscFunctionReturn(0); 2484 } 2485 2486 static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous) 2487 { 2488 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 2489 2490 PetscFunctionBegin; 2491 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2492 PetscValidPointer(continuous, 2); 2493 *continuous = lag->continuous; 2494 PetscFunctionReturn(0); 2495 } 2496 2497 static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous) 2498 { 2499 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data; 2500 2501 PetscFunctionBegin; 2502 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2503 lag->continuous = continuous; 2504 PetscFunctionReturn(0); 2505 } 2506 2507 /*@ 2508 PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity 2509 2510 Not Collective 2511 2512 Input Parameter: 2513 . sp - the PetscDualSpace 2514 2515 Output Parameter: 2516 . continuous - flag for element continuity 2517 2518 Level: intermediate 2519 2520 .seealso: PetscDualSpaceLagrangeSetContinuity() 2521 @*/ 2522 PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous) 2523 { 2524 PetscErrorCode ierr; 2525 2526 PetscFunctionBegin; 2527 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2528 PetscValidPointer(continuous, 2); 2529 ierr = PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));CHKERRQ(ierr); 2530 PetscFunctionReturn(0); 2531 } 2532 2533 /*@ 2534 PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous 2535 2536 Logically Collective on sp 2537 2538 Input Parameters: 2539 + sp - the PetscDualSpace 2540 - continuous - flag for element continuity 2541 2542 Options Database: 2543 . -petscdualspace_lagrange_continuity <bool> 2544 2545 Level: intermediate 2546 2547 .seealso: PetscDualSpaceLagrangeGetContinuity() 2548 @*/ 2549 PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous) 2550 { 2551 PetscErrorCode ierr; 2552 2553 PetscFunctionBegin; 2554 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2555 PetscValidLogicalCollectiveBool(sp, continuous, 2); 2556 ierr = PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));CHKERRQ(ierr); 2557 PetscFunctionReturn(0); 2558 } 2559 2560 static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Lagrange(PetscDualSpace sp, PetscBool *tensor) 2561 { 2562 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 2563 2564 PetscFunctionBegin; 2565 *tensor = lag->tensorSpace; 2566 PetscFunctionReturn(0); 2567 } 2568 2569 static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Lagrange(PetscDualSpace sp, PetscBool tensor) 2570 { 2571 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 2572 2573 PetscFunctionBegin; 2574 lag->tensorSpace = tensor; 2575 PetscFunctionReturn(0); 2576 } 2577 2578 static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Lagrange(PetscDualSpace sp, PetscBool *trimmed) 2579 { 2580 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 2581 2582 PetscFunctionBegin; 2583 *trimmed = lag->trimmed; 2584 PetscFunctionReturn(0); 2585 } 2586 2587 static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Lagrange(PetscDualSpace sp, PetscBool trimmed) 2588 { 2589 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 2590 2591 PetscFunctionBegin; 2592 lag->trimmed = trimmed; 2593 PetscFunctionReturn(0); 2594 } 2595 2596 static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent) 2597 { 2598 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 2599 2600 PetscFunctionBegin; 2601 if (nodeType) *nodeType = lag->nodeType; 2602 if (boundary) *boundary = lag->endNodes; 2603 if (exponent) *exponent = lag->nodeExponent; 2604 PetscFunctionReturn(0); 2605 } 2606 2607 static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Lagrange(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent) 2608 { 2609 PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *)sp->data; 2610 2611 PetscFunctionBegin; 2612 if (nodeType == PETSCDTNODES_GAUSSJACOBI && exponent <= -1.) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Exponent must be > -1"); 2613 lag->nodeType = nodeType; 2614 lag->endNodes = boundary; 2615 lag->nodeExponent = exponent; 2616 PetscFunctionReturn(0); 2617 } 2618 2619 /*@ 2620 PetscDualSpaceLagrangeGetTensor - Get the tensor nature of the dual space 2621 2622 Not collective 2623 2624 Input Parameter: 2625 . sp - The PetscDualSpace 2626 2627 Output Parameter: 2628 . tensor - Whether the dual space has tensor layout (vs. simplicial) 2629 2630 Level: intermediate 2631 2632 .seealso: PetscDualSpaceLagrangeSetTensor(), PetscDualSpaceCreate() 2633 @*/ 2634 PetscErrorCode PetscDualSpaceLagrangeGetTensor(PetscDualSpace sp, PetscBool *tensor) 2635 { 2636 PetscErrorCode ierr; 2637 2638 PetscFunctionBegin; 2639 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2640 PetscValidPointer(tensor, 2); 2641 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTensor_C",(PetscDualSpace,PetscBool *),(sp,tensor));CHKERRQ(ierr); 2642 PetscFunctionReturn(0); 2643 } 2644 2645 /*@ 2646 PetscDualSpaceLagrangeSetTensor - Set the tensor nature of the dual space 2647 2648 Not collective 2649 2650 Input Parameters: 2651 + sp - The PetscDualSpace 2652 - tensor - Whether the dual space has tensor layout (vs. simplicial) 2653 2654 Level: intermediate 2655 2656 .seealso: PetscDualSpaceLagrangeGetTensor(), PetscDualSpaceCreate() 2657 @*/ 2658 PetscErrorCode PetscDualSpaceLagrangeSetTensor(PetscDualSpace sp, PetscBool tensor) 2659 { 2660 PetscErrorCode ierr; 2661 2662 PetscFunctionBegin; 2663 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2664 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTensor_C",(PetscDualSpace,PetscBool),(sp,tensor));CHKERRQ(ierr); 2665 PetscFunctionReturn(0); 2666 } 2667 2668 /*@ 2669 PetscDualSpaceLagrangeGetTrimmed - Get the trimmed nature of the dual space 2670 2671 Not collective 2672 2673 Input Parameter: 2674 . sp - The PetscDualSpace 2675 2676 Output Parameter: 2677 . trimmed - Whether the dual space represents to dual basis of a trimmed polynomial space (e.g. Raviart-Thomas and higher order / other form degree variants) 2678 2679 Level: intermediate 2680 2681 .seealso: PetscDualSpaceLagrangeSetTrimmed(), PetscDualSpaceCreate() 2682 @*/ 2683 PetscErrorCode PetscDualSpaceLagrangeGetTrimmed(PetscDualSpace sp, PetscBool *trimmed) 2684 { 2685 PetscErrorCode ierr; 2686 2687 PetscFunctionBegin; 2688 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2689 PetscValidPointer(trimmed, 2); 2690 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetTrimmed_C",(PetscDualSpace,PetscBool *),(sp,trimmed));CHKERRQ(ierr); 2691 PetscFunctionReturn(0); 2692 } 2693 2694 /*@ 2695 PetscDualSpaceLagrangeSetTrimmed - Set the trimmed nature of the dual space 2696 2697 Not collective 2698 2699 Input Parameters: 2700 + sp - The PetscDualSpace 2701 - trimmed - Whether the dual space represents to dual basis of a trimmed polynomial space (e.g. Raviart-Thomas and higher order / other form degree variants) 2702 2703 Level: intermediate 2704 2705 .seealso: PetscDualSpaceLagrangeGetTrimmed(), PetscDualSpaceCreate() 2706 @*/ 2707 PetscErrorCode PetscDualSpaceLagrangeSetTrimmed(PetscDualSpace sp, PetscBool trimmed) 2708 { 2709 PetscErrorCode ierr; 2710 2711 PetscFunctionBegin; 2712 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2713 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetTrimmed_C",(PetscDualSpace,PetscBool),(sp,trimmed));CHKERRQ(ierr); 2714 PetscFunctionReturn(0); 2715 } 2716 2717 /*@ 2718 PetscDualSpaceLagrangeGetNodeType - Get a description of how nodes are laid out for Lagrange polynomials in this 2719 dual space 2720 2721 Not collective 2722 2723 Input Parameter: 2724 . sp - The PetscDualSpace 2725 2726 Output Parameters: 2727 + nodeType - The type of nodes 2728 . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that 2729 include the boundary are Gauss-Lobatto-Jacobi nodes) 2730 - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function 2731 '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type 2732 2733 Level: advanced 2734 2735 .seealso: PetscDTNodeType, PetscDualSpaceLagrangeSetNodeType() 2736 @*/ 2737 PetscErrorCode PetscDualSpaceLagrangeGetNodeType(PetscDualSpace sp, PetscDTNodeType *nodeType, PetscBool *boundary, PetscReal *exponent) 2738 { 2739 PetscErrorCode ierr; 2740 2741 PetscFunctionBegin; 2742 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2743 if (nodeType) PetscValidPointer(nodeType, 2); 2744 if (boundary) PetscValidPointer(boundary, 3); 2745 if (exponent) PetscValidPointer(exponent, 4); 2746 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeGetNodeType_C",(PetscDualSpace,PetscDTNodeType *,PetscBool *,PetscReal *),(sp,nodeType,boundary,exponent));CHKERRQ(ierr); 2747 PetscFunctionReturn(0); 2748 } 2749 2750 /*@ 2751 PetscDualSpaceLagrangeSetNodeType - Set a description of how nodes are laid out for Lagrange polynomials in this 2752 dual space 2753 2754 Logically collective 2755 2756 Input Parameters: 2757 + sp - The PetscDualSpace 2758 . nodeType - The type of nodes 2759 . boundary - Whether the node type is one that includes endpoints (if nodeType is PETSCDTNODES_GAUSSJACOBI, nodes that 2760 include the boundary are Gauss-Lobatto-Jacobi nodes) 2761 - exponent - If nodeType is PETSCDTNODES_GAUSJACOBI, indicates the exponent used for both ends of the 1D Jacobi weight function 2762 '0' is Gauss-Legendre, '-0.5' is Gauss-Chebyshev of the first type, '0.5' is Gauss-Chebyshev of the second type 2763 2764 Level: advanced 2765 2766 .seealso: PetscDTNodeType, PetscDualSpaceLagrangeGetNodeType() 2767 @*/ 2768 PetscErrorCode PetscDualSpaceLagrangeSetNodeType(PetscDualSpace sp, PetscDTNodeType nodeType, PetscBool boundary, PetscReal exponent) 2769 { 2770 PetscErrorCode ierr; 2771 2772 PetscFunctionBegin; 2773 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2774 ierr = PetscTryMethod(sp,"PetscDualSpaceLagrangeSetNodeType_C",(PetscDualSpace,PetscDTNodeType,PetscBool,PetscReal),(sp,nodeType,boundary,exponent));CHKERRQ(ierr); 2775 PetscFunctionReturn(0); 2776 } 2777 2778 2779 static PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp) 2780 { 2781 PetscFunctionBegin; 2782 sp->ops->destroy = PetscDualSpaceDestroy_Lagrange; 2783 sp->ops->view = PetscDualSpaceView_Lagrange; 2784 sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange; 2785 sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange; 2786 sp->ops->setup = PetscDualSpaceSetUp_Lagrange; 2787 sp->ops->createheightsubspace = NULL; 2788 sp->ops->createpointsubspace = NULL; 2789 sp->ops->getsymmetries = PetscDualSpaceGetSymmetries_Lagrange; 2790 sp->ops->apply = PetscDualSpaceApplyDefault; 2791 sp->ops->applyall = PetscDualSpaceApplyAllDefault; 2792 sp->ops->applyint = PetscDualSpaceApplyInteriorDefault; 2793 sp->ops->createalldata = PetscDualSpaceCreateAllDataDefault; 2794 sp->ops->createintdata = PetscDualSpaceCreateInteriorDataDefault; 2795 PetscFunctionReturn(0); 2796 } 2797 2798 /*MC 2799 PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals 2800 2801 Level: intermediate 2802 2803 .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType() 2804 M*/ 2805 PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp) 2806 { 2807 PetscDualSpace_Lag *lag; 2808 PetscErrorCode ierr; 2809 2810 PetscFunctionBegin; 2811 PetscValidHeaderSpecific(sp, PETSCDUALSPACE_CLASSID, 1); 2812 ierr = PetscNewLog(sp,&lag);CHKERRQ(ierr); 2813 sp->data = lag; 2814 2815 lag->tensorCell = PETSC_FALSE; 2816 lag->tensorSpace = PETSC_FALSE; 2817 lag->continuous = PETSC_TRUE; 2818 lag->numCopies = PETSC_DEFAULT; 2819 lag->numNodeSkip = PETSC_DEFAULT; 2820 lag->nodeType = PETSCDTNODES_DEFAULT; 2821 2822 ierr = PetscDualSpaceInitialize_Lagrange(sp);CHKERRQ(ierr); 2823 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);CHKERRQ(ierr); 2824 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);CHKERRQ(ierr); 2825 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Lagrange);CHKERRQ(ierr); 2826 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Lagrange);CHKERRQ(ierr); 2827 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Lagrange);CHKERRQ(ierr); 2828 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Lagrange);CHKERRQ(ierr); 2829 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Lagrange);CHKERRQ(ierr); 2830 ierr = PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Lagrange);CHKERRQ(ierr); 2831 PetscFunctionReturn(0); 2832 } 2833 2834