1 #include <petsc/private/dmfieldimpl.h> /*I "petscdmfield.h" I*/ 2 #include <petscfe.h> 3 #include <petscdmplex.h> 4 #include <petscds.h> 5 6 typedef struct _n_DMField_DS 7 { 8 PetscInt fieldNum; 9 Vec vec; 10 PetscInt height; 11 PetscObject *disc; 12 PetscBool multifieldVec; 13 } 14 DMField_DS; 15 16 static PetscErrorCode DMFieldDestroy_DS(DMField field) 17 { 18 DMField_DS *dsfield; 19 PetscInt i; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 dsfield = (DMField_DS *) field->data; 24 ierr = VecDestroy(&dsfield->vec);CHKERRQ(ierr); 25 for (i = 0; i < dsfield->height; i++) { 26 ierr = PetscObjectDereference(dsfield->disc[i]);CHKERRQ(ierr); 27 } 28 ierr = PetscFree(dsfield->disc);CHKERRQ(ierr); 29 ierr = PetscFree(dsfield);CHKERRQ(ierr); 30 PetscFunctionReturn(0); 31 } 32 33 static PetscErrorCode DMFieldView_DS(DMField field,PetscViewer viewer) 34 { 35 DMField_DS *dsfield = (DMField_DS *) field->data; 36 PetscBool iascii; 37 PetscObject disc; 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 42 disc = dsfield->disc[0]; 43 if (iascii) { 44 PetscViewerASCIIPrintf(viewer, "PetscDS field %D\n",dsfield->fieldNum);CHKERRQ(ierr); 45 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 46 ierr = PetscObjectView(disc,viewer);CHKERRQ(ierr); 47 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 48 } 49 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 50 if (dsfield->multifieldVec) { 51 SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"View of subfield not implemented yet"); 52 } else { 53 ierr = VecView(dsfield->vec,viewer);CHKERRQ(ierr); 54 } 55 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 56 PetscFunctionReturn(0); 57 } 58 59 static PetscErrorCode DMFieldEvaluate_DS(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H) 60 { 61 PetscFunctionBegin; 62 SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented yet"); 63 PetscFunctionReturn(0); 64 } 65 66 static PetscErrorCode DMFieldDSGetHeightDisc(DMField field, PetscInt height, PetscObject *disc) 67 { 68 DMField_DS *dsfield = (DMField_DS *) field->data; 69 PetscErrorCode ierr; 70 71 PetscFunctionBegin; 72 if (!dsfield->disc[height]) { 73 PetscClassId id; 74 75 ierr = PetscObjectGetClassId(dsfield->disc[0],&id);CHKERRQ(ierr); 76 if (id == PETSCFE_CLASSID) { 77 PetscFE fe = (PetscFE) dsfield->disc[0]; 78 79 ierr = PetscFECreateHeightTrace(fe,height,(PetscFE *)&dsfield->disc[height]);CHKERRQ(ierr); 80 } 81 } 82 *disc = dsfield->disc[height]; 83 PetscFunctionReturn(0); 84 } 85 86 #define DMFieldDSdot(y,A,b,m,n,c,cast) \ 87 do { \ 88 PetscInt _i, _j, _k; \ 89 for (_i = 0; _i < (m); _i++) { \ 90 for (_k = 0; _k < (c); _k++) { \ 91 (y)[_i * (c) + _k] = 0.; \ 92 } \ 93 for (_j = 0; _j < (n); _j++) { \ 94 for (_k = 0; _k < (c); _k++) { \ 95 (y)[_i * (c) + _k] += (A)[(_i * (n) + _j) * (c) + _k] * cast((b)[_j]); \ 96 } \ 97 } \ 98 } \ 99 } while (0) 100 101 static PetscErrorCode DMFieldEvaluateFE_DS(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H) 102 { 103 DMField_DS *dsfield = (DMField_DS *) field->data; 104 DM dm; 105 PetscObject disc; 106 PetscClassId classid; 107 PetscInt nq, nc, dim, meshDim, numCells; 108 PetscSection section; 109 const PetscReal *qpoints; 110 PetscBool isStride; 111 const PetscInt *points = NULL; 112 PetscInt sfirst = -1, stride = -1; 113 PetscErrorCode ierr; 114 115 PetscFunctionBeginHot; 116 dm = field->dm; 117 nc = field->numComponents; 118 ierr = PetscQuadratureGetData(quad,&dim,NULL,&nq,&qpoints,NULL);CHKERRQ(ierr); 119 ierr = DMFieldDSGetHeightDisc(field,dsfield->height - dim,&disc);CHKERRQ(ierr); 120 ierr = DMGetDimension(dm,&meshDim);CHKERRQ(ierr); 121 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 122 ierr = PetscSectionGetField(section,dsfield->fieldNum,§ion);CHKERRQ(ierr); 123 ierr = PetscObjectGetClassId(disc,&classid);CHKERRQ(ierr); 124 /* TODO: batch */ 125 ierr = PetscObjectTypeCompare((PetscObject)pointIS,ISSTRIDE,&isStride);CHKERRQ(ierr); 126 ierr = ISGetLocalSize(pointIS,&numCells);CHKERRQ(ierr); 127 if (isStride) { 128 ierr = ISStrideGetInfo(pointIS,&sfirst,&stride);CHKERRQ(ierr); 129 } else { 130 ierr = ISGetIndices(pointIS,&points);CHKERRQ(ierr); 131 } 132 if (classid == PETSCFE_CLASSID) { 133 PetscFE fe = (PetscFE) disc; 134 PetscInt feDim, i; 135 PetscReal *fB = NULL, *fD = NULL, *fH = NULL; 136 PetscInt closureSize; 137 PetscScalar *elem = NULL; 138 139 if (dim == meshDim - 1) { 140 /* TODO */ 141 } 142 ierr = PetscFEGetDimension(fe,&feDim);CHKERRQ(ierr); 143 ierr = PetscFEGetTabulation(fe,nq,qpoints,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);CHKERRQ(ierr); 144 closureSize = feDim; 145 for (i = 0; i < numCells; i++) { 146 PetscInt c = isStride ? (sfirst + i * stride) : points[i]; 147 148 ierr = DMPlexVecGetClosure(dm,section,dsfield->vec,c,&closureSize,&elem);CHKERRQ(ierr); 149 if (B) { 150 if (type == PETSC_SCALAR) { 151 PetscScalar *cB = &((PetscScalar *) B)[nc * nq * i]; 152 153 DMFieldDSdot(cB,fB,elem,nq,feDim,nc,(PetscScalar)); 154 } else { 155 PetscReal *cB = &((PetscReal *) B)[nc * nq * i]; 156 157 DMFieldDSdot(cB,fB,elem,nq,feDim,nc,PetscRealPart); 158 } 159 } 160 if (D) { 161 if (type == PETSC_SCALAR) { 162 PetscScalar *cD = &((PetscScalar *) D)[nc * nq * dim * i]; 163 164 DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),(PetscScalar)); 165 } else { 166 PetscReal *cD = &((PetscReal *) D)[nc * nq * dim * i]; 167 168 DMFieldDSdot(cD,fD,elem,nq,feDim,(nc * dim),PetscRealPart); 169 } 170 } 171 if (H) { 172 if (type == PETSC_SCALAR) { 173 PetscScalar *cH = &((PetscScalar *) H)[nc * nq * dim * dim * i]; 174 175 DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),(PetscScalar)); 176 } else { 177 PetscReal *cH = &((PetscReal *) H)[nc * nq * dim * dim * i]; 178 179 DMFieldDSdot(cH,fH,elem,nq,feDim,(nc * dim * dim),PetscRealPart); 180 } 181 } 182 } 183 ierr = DMRestoreWorkArray(dm,feDim,MPIU_SCALAR,&elem);CHKERRQ(ierr); 184 ierr = PetscFERestoreTabulation(fe,nq,qpoints,B ? &fB : NULL,D ? &fD : NULL,H ? &fH : NULL);CHKERRQ(ierr); 185 } else {SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented");} 186 if (!isStride) { 187 ierr = ISRestoreIndices(pointIS,&points);CHKERRQ(ierr); 188 } 189 ierr = PetscSectionDestroy(§ion);CHKERRQ(ierr); 190 PetscFunctionReturn(0); 191 } 192 193 static PetscErrorCode DMFieldGetFEInvariance_DS(DMField field, IS pointIS, PetscBool *isConstant, PetscBool *isAffine, PetscBool *isQuadratic) 194 { 195 DMField_DS *dsfield; 196 PetscObject disc; 197 PetscInt h, imin; 198 PetscClassId id; 199 PetscErrorCode ierr; 200 201 PetscFunctionBegin; 202 dsfield = (DMField_DS *) field->data; 203 ierr = ISGetMinMax(pointIS,&imin,NULL);CHKERRQ(ierr); 204 for (h = 0; h <= dsfield->height; h++) { 205 PetscInt hEnd; 206 207 ierr = DMPlexGetHeightStratum(field->dm,h,NULL,&hEnd);CHKERRQ(ierr); 208 if (imin < hEnd) break; 209 } 210 ierr = DMFieldDSGetHeightDisc(field,h,&disc);CHKERRQ(ierr); 211 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 212 if (id == PETSCFE_CLASSID) { 213 PetscFE fe = (PetscFE) disc; 214 PetscInt order, maxOrder; 215 PetscBool tensor = PETSC_FALSE; 216 PetscSpace sp; 217 218 ierr = PetscFEGetBasisSpace(fe, &sp);CHKERRQ(ierr); 219 ierr = PetscSpaceGetOrder(sp,&order);CHKERRQ(ierr); 220 ierr = PetscSpacePolynomialGetTensor(sp,&tensor);CHKERRQ(ierr); 221 if (tensor) { 222 PetscInt dim; 223 224 ierr = DMGetDimension(field->dm,&dim);CHKERRQ(ierr); 225 maxOrder = order * dim; 226 } else { 227 maxOrder = order; 228 } 229 if (isConstant) *isConstant = (maxOrder < 1) ? PETSC_TRUE : PETSC_FALSE; 230 if (isAffine) *isAffine = (maxOrder < 2) ? PETSC_TRUE : PETSC_FALSE; 231 if (isQuadratic) *isQuadratic = (maxOrder < 3) ? PETSC_TRUE : PETSC_FALSE; 232 } 233 PetscFunctionReturn(0); 234 } 235 236 static PetscErrorCode DMFieldCreateDefaultQuadrature_DS(DMField field, IS pointIS, PetscQuadrature *quad) 237 { 238 PetscInt h, dim, imax, imin; 239 DM dm; 240 DMField_DS *dsfield; 241 PetscObject disc; 242 PetscFE fe; 243 PetscClassId id; 244 PetscErrorCode ierr; 245 246 247 PetscFunctionBegin; 248 dm = field->dm; 249 dsfield = (DMField_DS *) field->data; 250 ierr = ISGetMinMax(pointIS,&imax,&imin);CHKERRQ(ierr); 251 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr); 252 for (h = 0; h <= dim; h++) { 253 PetscInt hStart, hEnd; 254 255 ierr = DMPlexGetHeightStratum(dm,h,&hStart,&hEnd);CHKERRQ(ierr); 256 if (imin >= hStart && imax < hEnd) break; 257 } 258 *quad = NULL; 259 if (h <= dsfield->height) { 260 ierr = DMFieldDSGetHeightDisc(field,h,&disc);CHKERRQ(ierr); 261 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 262 if (id != PETSCFE_CLASSID) PetscFunctionReturn(0); 263 fe = (PetscFE) disc; 264 ierr = PetscFEGetQuadrature(fe,quad);CHKERRQ(ierr); 265 ierr = PetscObjectReference((PetscObject)*quad);CHKERRQ(ierr); 266 } 267 PetscFunctionReturn(0); 268 } 269 270 static PetscErrorCode DMFieldInitialize_DS(DMField field) 271 { 272 PetscFunctionBegin; 273 field->ops->destroy = DMFieldDestroy_DS; 274 field->ops->evaluate = DMFieldEvaluate_DS; 275 field->ops->evaluateFE = DMFieldEvaluateFE_DS; 276 field->ops->getFEInvariance = DMFieldGetFEInvariance_DS; 277 field->ops->createDefaultQuadrature = DMFieldCreateDefaultQuadrature_DS; 278 field->ops->view = DMFieldView_DS; 279 PetscFunctionReturn(0); 280 } 281 282 PETSC_INTERN PetscErrorCode DMFieldCreate_DS(DMField field) 283 { 284 DMField_DS *dsfield; 285 PetscErrorCode ierr; 286 287 PetscFunctionBegin; 288 ierr = PetscNewLog(field,&dsfield);CHKERRQ(ierr); 289 field->data = dsfield; 290 ierr = DMFieldInitialize_DS(field);CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 PetscErrorCode DMFieldCreateDS(DM dm, PetscInt fieldNum, Vec vec,DMField *field) 295 { 296 DMField b; 297 DMField_DS *dsfield; 298 PetscObject disc = NULL; 299 PetscBool isContainer = PETSC_FALSE; 300 PetscClassId id = -1; 301 PetscInt numComponents = -1, dsNumFields; 302 PetscSection section; 303 PetscErrorCode ierr; 304 305 PetscFunctionBegin; 306 ierr = DMGetDefaultSection(dm,§ion);CHKERRQ(ierr); 307 ierr = PetscSectionGetFieldComponents(section,fieldNum,&numComponents);CHKERRQ(ierr); 308 ierr = DMGetNumFields(dm,&dsNumFields);CHKERRQ(ierr); 309 if (dsNumFields) {ierr = DMGetField(dm,fieldNum,&disc);CHKERRQ(ierr);} 310 if (disc) { 311 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 312 isContainer = (id == PETSC_CONTAINER_CLASSID) ? PETSC_TRUE : PETSC_FALSE; 313 } 314 if (!disc || isContainer) { 315 MPI_Comm comm = PetscObjectComm((PetscObject) dm); 316 PetscInt cStart, cEnd, dim; 317 PetscInt localConeSize = 0, coneSize; 318 PetscFE fe; 319 PetscDualSpace Q; 320 PetscSpace P; 321 DM K; 322 PetscQuadrature quad, fquad; 323 PetscBool isSimplex; 324 325 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 326 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 327 if (cEnd > cStart) { 328 ierr = DMPlexGetConeSize(dm, cStart, &localConeSize);CHKERRQ(ierr); 329 } 330 ierr = MPI_Allreduce(&localConeSize,&coneSize,1,MPIU_INT,MPI_MAX,comm);CHKERRQ(ierr); 331 isSimplex = (coneSize == (dim + 1)) ? PETSC_TRUE : PETSC_FALSE; 332 ierr = PetscSpaceCreate(comm, &P);CHKERRQ(ierr); 333 ierr = PetscSpaceSetOrder(P, 1);CHKERRQ(ierr); 334 ierr = PetscSpaceSetNumComponents(P, numComponents);CHKERRQ(ierr); 335 ierr = PetscSpaceSetType(P,PETSCSPACEPOLYNOMIAL);CHKERRQ(ierr); 336 ierr = PetscSpaceSetNumVariables(P, dim);CHKERRQ(ierr); 337 ierr = PetscSpacePolynomialSetTensor(P, isSimplex ? PETSC_FALSE : PETSC_TRUE);CHKERRQ(ierr); 338 ierr = PetscSpaceSetUp(P);CHKERRQ(ierr); 339 ierr = PetscDualSpaceCreate(comm, &Q);CHKERRQ(ierr); 340 ierr = PetscDualSpaceSetType(Q,PETSCDUALSPACELAGRANGE);CHKERRQ(ierr); 341 ierr = PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);CHKERRQ(ierr); 342 ierr = PetscDualSpaceSetDM(Q, K);CHKERRQ(ierr); 343 ierr = DMDestroy(&K);CHKERRQ(ierr); 344 ierr = PetscDualSpaceSetNumComponents(Q, numComponents);CHKERRQ(ierr); 345 ierr = PetscDualSpaceSetOrder(Q, 1);CHKERRQ(ierr); 346 ierr = PetscDualSpaceLagrangeSetTensor(Q, isSimplex ? PETSC_FALSE : PETSC_TRUE);CHKERRQ(ierr); 347 ierr = PetscDualSpaceSetUp(Q);CHKERRQ(ierr); 348 ierr = PetscFECreate(comm, &fe);CHKERRQ(ierr); 349 ierr = PetscFESetType(fe,PETSCFEBASIC);CHKERRQ(ierr); 350 ierr = PetscFESetBasisSpace(fe, P);CHKERRQ(ierr); 351 ierr = PetscFESetDualSpace(fe, Q);CHKERRQ(ierr); 352 ierr = PetscFESetNumComponents(fe, numComponents);CHKERRQ(ierr); 353 ierr = PetscFESetUp(fe);CHKERRQ(ierr); 354 ierr = PetscSpaceDestroy(&P);CHKERRQ(ierr); 355 ierr = PetscDualSpaceDestroy(&Q);CHKERRQ(ierr); 356 if (isSimplex) { 357 ierr = PetscDTGaussJacobiQuadrature(dim, 1, 1, -1.0, 1.0, &quad);CHKERRQ(ierr); 358 ierr = PetscDTGaussJacobiQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);CHKERRQ(ierr); 359 } 360 else { 361 ierr = PetscDTGaussTensorQuadrature(dim, 1, 1, -1.0, 1.0, &quad);CHKERRQ(ierr); 362 ierr = PetscDTGaussTensorQuadrature(dim-1, 1, 1, -1.0, 1.0, &fquad);CHKERRQ(ierr); 363 } 364 ierr = PetscFESetQuadrature(fe, quad);CHKERRQ(ierr); 365 ierr = PetscFESetFaceQuadrature(fe, fquad);CHKERRQ(ierr); 366 ierr = PetscQuadratureDestroy(&quad);CHKERRQ(ierr); 367 ierr = PetscQuadratureDestroy(&fquad);CHKERRQ(ierr); 368 disc = (PetscObject) fe; 369 } else { 370 ierr = PetscObjectReference(disc);CHKERRQ(ierr); 371 } 372 ierr = PetscObjectGetClassId(disc,&id);CHKERRQ(ierr); 373 if (id == PETSCFE_CLASSID) { 374 PetscFE fe = (PetscFE) disc; 375 376 ierr = PetscFEGetNumComponents(fe,&numComponents);CHKERRQ(ierr); 377 } else {SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented");} 378 ierr = DMFieldCreate(dm,numComponents,DMFIELD_VERTEX,&b);CHKERRQ(ierr); 379 ierr = DMFieldSetType(b,DMFIELDDS);CHKERRQ(ierr); 380 dsfield = (DMField_DS *) b->data; 381 dsfield->fieldNum = fieldNum; 382 ierr = DMGetDimension(dm,&dsfield->height);CHKERRQ(ierr); 383 ierr = PetscCalloc1(dsfield->height,&dsfield->disc);CHKERRQ(ierr); 384 dsfield->disc[0] = disc; 385 ierr = PetscObjectReference((PetscObject)vec);CHKERRQ(ierr); 386 dsfield->vec = vec; 387 *field = b; 388 389 PetscFunctionReturn(0); 390 } 391